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uncertainties  would  cause  significantly  different  results.  A  significant  shift  to  larger  bubbles  and 
larger  swelling  was  indicated  near  1800°K. 

Poorly  defined  property  values  were  varied  and  results  compared.  The  largest  effects 
were  observed  from  the  irradiation  enhanced  diffusion  coefficient  and  the  resolution  parameter. 
The  resolution  parameter  uncertainty  had  the  largest  potential  effect  on  results. 
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ABSTRACT 

Transient  Fission  Gas  Behavior  in  Uranium  Nitride  Fuel 
Under  Proposed  Space  Applications.  (December  1991) 

Daniel  Lee  DeForest,  B.S.,  Baylor  University; 

M  S.,  Air  Force  Institute  of  Technology 
Chair  of  Advisory  Committee;  Dr.  K.L.  Peddicord 

In  order  to  investigate  whether  fission  gas  swelling  and  release  would  be 
significant  factors  in  a  space  based  nuclear  reactor  operating  under  the  Strategic 
Defense  Initiative  (SDI)  program,  the  finite  element  program  REDSTONE  (Routine 
for  Evaluating  Dynamic  ^welling  in  Transient  Operational  Nuclear  Environments) 
was  developed  to  model  the  1-D,  spherical  geometry  diffusion  equations  describing 
transient  fission  gas  behavior  in  a  single  uranium  nitride  fuel  grain.  The  equations 
characterized  individual  bubbles,  rather  than  bubble  groupings.  This  limited 
calculations  to  those  scenarios  where  low  temperatures,  low  bumups,  or  both  were 
present.  Instabilities  in  the  bubble  radii  calculations  forced  the  implementation  of 
additional  constraints  limiting  the  bubble  sizes  to  minimum  and  maximum 
(equilibrium)  radii. 

The  validity  of  REDSTONE  calculations  were  checked  against  analytical 
solutions  for  internal  consistency  and  agaiinst  experimental  studies  for  agreement 
with  swelling  and  release  results.  These  checks  indicated  that  REDSTONE  was 
satisfactorily  accounting  for  the  intragranular  fission  gas  portion  of  swelling,  but  that 
other  non-gaseous  sources  accounted  for  up  to  1.5%  additional  swelling. 


IV 


A  calculation  was  performed  to  determine  the  swelling  and  gas  release 
associated  with  a  hypothetical  SDI  scenario.  The  SDI  scenario  involved  10  years  of 
low  power  operation  at  10  kW  followed  by  a  30  second  rise  to  1  MW  which  was 
maintained  for  approximately  15  minutes.  The  results  indicated  that  the  fuel  pellet 
diametral  increase  would  remain  less  than  1%  and  release  would  be  negligible,  far 
below  that  needed  for  design  failure.  This  same  scenario  was  run  at  higher 
temperatures  to  see  if  uncertainties  in  the  operating  parameters  could  cause 
significantly  different  results.  It  was  apparent  that  further  investigations  into  the 
swelling  behavior  near  1 800“K  were  needed,  as  results  indicated  a  significant  shift 
to  larger  bubbles  at  this  temperature. 

In  an  effort  to  understand  the  effect  of  property  uncertainties  on  the  solution, 
several  of  the  properties  were  varied  and  the  results  compared.  The  largest  effects 
were  observed  from  changes  in  the  irradiation  enhanced  diffusion  coefficient  and  the 
resolution  parameter.  The  resolution  parameter  uncertainty  would  have  the  largest 
potential  effect  on  results,  likely  resulting  in  a  reduction  in  the  swelling  and  a 
corresponding  increase  in  release. 
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CHAPTER  I 
INTRODUCTION 


1.1  Justification 

Proposed  devices  for  use  in  the  Strategic  Defense  Initiative  (SDI)  program 
will  require  at  least  a  megawatt  of  power  and  will  operate  under  stringent 
requirements.  In  addition  to  the  harsh  environmental  requirements  of  space,  the 
operational  requirements  will  be  quite  severe.  These  latter  requirements  are  the 
ability  to  survive  in  a  low  power  ‘housekeeping’  mode  for  an  extended  period, 
perhaps  up  to  ten  years,  with  the  ability  to  come  to  full  power  operation  within 
seconds  and  maintain  this  power  for  a  minimum  of  several  minutes.  Nuclear  power 
is  the  most  viable  power  option  because  of  the  constraints  of  power/weight  ratio  and 
hands-off  operation. 

Specific  operational  parameters  for  a  multi-megawatt,  SDI  reactor  are  not 
publicly  available.  However,  one  can  imagine  a  piossible  scenario  where  a  reactor  is 
sitting  in  ‘housekeeping’  mode  at  just  enough  power  to  keep  the  system  functional, 
perhaps  10  kW  electric.  Upon  command,  a  rapid  increase  to  full  power  of  1  MW 
electric  would  then  be  executed,  where  the  system  would  stay  for  several  minutes. 
Conceivably,  the  ramp  to  full  power  might  occur  on  the  order  of  30  seconds  with  the 
full  power  mode  lasting  for  15  minutes.  The  full  power  increase  from  10  kW  to 
IMW  in  just  30  seconds  could  place  severe  thermal  transients  on  the  fuel.  The  high 
power  operation  could  cause  a  large  inventory  of  fission  gas  to  coalesce  into  large 
bubbles  inducing  unwanted  swelling  of  the  fuel,  or  it  could  result  in  quick  release  of 
a  large  quantity  of  fission  gas  producing  unwanted  pressures  in  the  fuel  pin. 


Style  and  format  of  this  dissertation  follow  the  Journal  of  Nuclear  Materials. 


One  of  the  primary  design  constraints  of  an  SDI  nuclear  system  will  be  the 
amount  of  pin  diameter  increase  allowed  before  the  system  is  determined  to  have 
‘failed’.  An  increase  in  fuel  pin  diameter  occurs  either  directly  from  expansion  of 
the  fuel  j^ainst  the  cladding,  or  through  excess  internal  pressure  from  fission  gas 
which  is  released  to  the  fuel-clad  gap.  The  former  case  depends  strongly  on  local 
conditions  of  the  fuel  pin  since  the  fission  gas  is  ‘trapped’  in  the  fuel  at  a  specific 
location.  In  this  case,  the  fission  gas  behavior  both  within  the  fuel  grains 
(intragranular  fission  gas)  and  on  the  grain  boundaries  (intergranular  fission  gas) 
must  be  determined.  The  latter  case  depends  less  on  local  conditions  since  a  large 
quantity  of  released  fission  gas  will  create  open  pathways  out  of  the  fuel  and  much 
of  the  gas  that  reaches  the  grain  boundaries  will  be  released  to  the  gap  where  it 
mixes  with  gas  released  from  other  fuel  pin  regions.  Here,  intragranular  fission  gas 
behavior  is  of  primary  importance  since  it  is  implied  by  release  that  intergranular 
gas  does  not  play  a  significant  role.  It  is  important  to  realize,  however,  that  the 
actual  behavior  of  the  fuel  will  be  some  combination  of  the  above  two  extremes. 
Nuclear  fuel  is  usually  designed  to  either  release  the  fission  gas  to  the  gap  or  to 
accommodate  the  fission  gas  in  the  fuel.  In  either  case,  it  is  clear  that  the 
intragranular  fission  gas  behavior  must  first  be  known  in  order  to  properly  design 
the  fuel.  The  importance  of  considering  intergranular  fission  gas  behavior  is 
determined  by  the  specific  design  of  the  fuel  used. 

The  effects  of  thermally  expanding  fuel  (fission  gas  swelling  excluded)  and 
the  associated  stresses  and  strains  can  be  determined  from  macroscopic  properties  of 
the  fuel  and  thus  are  easier  to  include  in  a  design.  However,  the  effects  from  fission 
gas,  requiring  knowledge  of  the  microscopic  properties  of  the  fuel,  are  harder  to 
determine.  For  the  current  scenario,  the  initial  conditions  of  long  term,  low  power 
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operation  followed  by  a  severe  increase  in  power  by  two  orders  of  magnitude  are 
unique  to  this  system.  A  prediction  of  possible  fission  gas  behavior  in  such  a  system 
has  not  been  previously  reported.  Potentially,  the  fission  gas  could  simply  be 
released  from  the  grain  to  grain  boundaries  which  would  become  open  avenues  to 
the  plenum  either  through  cracks  or  interconnected  grain  boundary  tunnels 
developed  during  the  transient.  On  the  other  hand,  the  fission  gas  could  be 
completely  retained  by  the  fuel,  i.e.  no  cracks,  and  expansion  of  retained  fission  gas 
as  bubbles  could  cause  severe  swelling  of  the  fuel.  The  severity  of  this  additional 
swelling  would  depend  on  where  the  fission  gas  resided  during  the  transient. 
Intragranular  fission  gas  would  be  present  either  as  individual  gas  atoms  or  in  vziy 
small  bubbles.  Intergranular  fission  gas  tends  to  form  much  larger,  lower  density 
bubbles  than  found  inside  the  grain.  Grain  boundaries  are  not  as  ‘stiff  as  the 
interior  of  the  grain  to  expansion  by  excess  bubble  pressure.  As  a  result,  fission  gas 
retained  by  the  fuel  on  the  grain  boundaries  will  cause  the  most  severe  swelling  of 
the  fuel.  When  intergranular  fission  gas  is  retained,  intragranular  fission  gas  plays  a 
lesser  role  in  direct  swelling  of  the  fuel.  However,  knowledge  of  intragranular 
fission  gas  behavior  is  important  as  it  is  a  feed  for  both  intergranular  bubbles  and 
gas  release.  Thus,  one  must  first  analyze  the  intragranular  fission  gas  behavior  as  a 
prerequisite  to  determining  the  effects  from  other  fission  gas  processes. 

1.2  Initial  Indications 

Scoping  calculations  based  on  the  system  and  scenario  as  outlined  above 
provided  some  insight  into  which  processes  would  need  to  be  included  in  modeling 
UN  fuel.  During  the  low  power  phase  of  the  scenario,  the  fuel  temperature  appears 
to  be  cool  enough  that  fission  gas  atoms  are  immobile  in  the  lattice.  Also,  there  is 
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only  a  small  temperature  gradient  in  the  fuel  due  to  the  higher  thermal  conductivity 
of  UN.  Bubbles  that  do  form  will  be  from  random  diffusion  with  additional  effects 
from  the  high  density  of  fission  gas  atoms  dejxjsited  over  a  ten  year  period.  It  is 
expected  that  during  this  phase  the  fission  gas  remains  in  small  bubbles  and  single 
atoms  and  within  the  fuel  grain.  During  the  high  power  phase  of  the  scenario,  the 
fuel  temperature  and  temperature  gradient  increases  so  that  gas  atoms  and  bubbles 
become  more  mobile.  However,  the  time  at  full  power  is  short  enough  that  any 
bubble  movement  is  severely  limited.  Initial  calculations  indicated  that  some  release 
from  the  fuel  grain  and  some  bubble  size  increase  might  occur,  but  not  to  an  extent 
that  would  cause  fuel  pin  failure.  In  summary,  initial  indications  were  that  the 
scenario  being  modeled  would  not  cause  large  swelling  or  large  release  from  the  fuel 
grain,  and  that  most  of  the  fission  gas  would  be  retained  in  the  grain.  However, 
uncertainties  in  experimentally  measured  fuel  properties  for  UN  will  cause 
corresponding  uncertainties  in  any  fission  gas  results  calculated. 

U  Use  of  Diffusion-Only  Equations 

As  a  result  of  scoping  calculations,  the  computer  program  REDSTONE 
(Routine  for  Evaluation  of  Dynamic  ^welling  in  Transient  Operational  Nuclear 
Environments),  was  written.  This  program  is  based  on  the  diffusion  equation  for  a 
1-D  spherical  fuel  grain.  Effects  from  thermal  gradients  are  included  as  additional 
coalescence  terms  to  the  diffusion  equation  and  in  some  cases  as  simple  release 
fractions  to  maintain  spherical  symmetry.  This  latter  approach  appears  justified  as 
the  expected  result  is  that  the  time  at  high  power  will  He  too  short  for  significant  gas 
movement  to  destroy  spherical  symmetry.  This  is  fortunate  since  detailed  modeling 
of  asymmetric  effects  would  cause  the  subsequent  profiles  within  the  fuel  grain  to 
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become  non-symmetric  and  necessitate  difficult  and  unwieldy  3-D  modeling  of  the 
fuel  grain. 

1.4  Use  of  Finite  Elements 

Program  REDSTONE  was  developed  as  a  finite  element  code  based  on  the 
1-D  non-linear  diffusion  equation  in  spherical  coordinates.  To  be  able  to  analyze  the 
non-linear  equation,  the  equations  had  to  be  incrementalized  and  solved  through 
iteration.  This  technique  has  been  used  previously  for  structural  mechanics. 
However,  the  type  of  non-linearity  present  in  the  diffusion  equation  used  here  is 
altogether  different  from  that  found  in  structural  mechanics.  In  structural  mechanics 
the  non-linearities  occur  both  as  non-constant  material  properties  and  as  cross 
products  in  the  displacement  gradients.  In  the  formulation  of  the  diffusion  equation 
used  here,  the  material  properties  are  held  constant  (i.e.  they  are  assumed  not  to  be  a 
fimction  of  fission  gas  atom  concentration)  and  the  non-linearities  occur  as  cross 
products  in  the  concentrations  themselves,  not  in  the  concentration  gradients.  Thus 
the  finite  element  formulation  of  the  diffusion  equation  which  was  used  in 
REDSTONE  is  fundamentally  different  from  that  used  in  structural  mechanics.  It 
appears  that  this  approach  has  not  been  reported  previously. 

1.5  Basic  Structure  of  REDSTONE 

The  decision  was  made  after  some  investigation  into  the  formulation  of  the 
equations  that  the  modeling  of  vacancies/interstitials  within  the  fuel  matrix  and 
bubble  radii  should  be  decoupled  from  the  calculation  of  bubble  concentrations.  It 
was  found  that  inclusion  of  the  vacancy/interstitial  equations  and  bubble  radii 
equations  with  concentration  equations  would  cause  such  a  high  degree  of 
non-linearity  that  an  accurate  finite  element  formulation  using  the  methods  applied 


6 


here  could  not  be  developed.  The  higher  non-linearity  would  have  reduced  the 
concentration  equations  down  to  an  essentially  linear  form  destroying  the  benefits  of 
using  the  method  in  the  first  place.  As  a  result,  the  sequence  of  calculations  in 
REDSTONE  is  to  first  calculate  concentration  changes,  keeping  vacancy/interstitial 
concentrations  and  bubble  radii  constant.  This  is  followed  by  calculation  of  bulk 
vacancy  and  interstitial  changes,  keeping  bubble  concentrations  and  radii  constant. 
After  calculation  of  vacancy/interstitial  concentrations,  radii  are  recomputed  and  the 
entire  process  repeated  until  no  additional  iterations  are  required  in  the  individual 
steps  of  bubble  concentration  calculation  and  radii  calculation.  Once  convergence 
as  described  above  has  occurred,  results  are  output,  the  time  is  incremented,  and  the 
process  starts  anew  until  the  end  of  the  problem.  The  calculation  of  bulk 
vacancv/interstitial  concentrations  and  bubble  radii  are  separated  as  the  latter 
calculation  is  strictly  linear  and  thus  more  efficient  to  calculate  alone. 

1.6  Objectives 

The  current  work  attempts  to  add  insight  to  the  modeling  of  fission  gas 
within  UN  and  answer  some  of  the  questions  posed  by  the  previously  described 
scenario.  The  objectives  of  this  study  are; 

a)  to  investigate  the  modeling  of  fission  gas  behavior  in  uranium  nitride 
(UN) 

b)  to  investigate  the  modeling  of  fission  gas  rate  equations  via  a  complete 
finite  element  formulation 

c)  the  determination  of  how  much  intragranular  fission  gas  migrates  to  the 
grain  boundaries  and  the  amount  of  intragranular  swelling 

d)  the  effect  of  uncertainties  in  the  currently  available  property  data 


7 


CHAPTER  II 

SCOPING  CALCULATIONS 

2.1  Choice  of  System  Parameters 

2.1.1  Selection  of  a  Nuclear  System 

Due  to  the  lack  of  a  current  established  design  for  an  SDI  nuclear  system  as 
well  as  the  classified  nature  of  the  SDI  program,  a  suitable  substitute  system  needed 
to  be  chosen  for  this  study.  The  most  obvious  substitute  was  the  SP- 1 00  nuclear 
reactor  proposed  for  electric  power  requirements  of  space  vehicles  and  satellites. 
This  reactor  system  has  an  established  design  and  it  is  conceivable  that  testing  of 
SDI  devices  would  utilize  the  SP-100  platform  first  [1].  Currently,  the  SP-100  is 
proposed  as  the  most  viable  power  option  for  systems  requiring  10  kW  to  1  MW  of 
electric  power.  The  figure-of- merit  for  failure  of  the  fuel  elements  in  this  reactor 
has  also  been  established.  This  limit  is  one  percent  diametral  increase  in  the  fuel 
element.  Essential  physical  characteristics  of  SP-100  in  the  1  MW  configuration  are 
listed  in  table  1. 

2.1.2  Selection  of  a  Nuclear  Fuel 

Because  of  the  long  lifetime  required  of  the  SP-100  reactor,  as  well  as  the 
need  to  have  a  compact,  lightweight,  and  non-refuelable  system,  UN  was  chosen  as 
the  nuclear  fuel  for  this  reactor.  The  primary  advantages  of  UN  over  UO2  and  UC 
include  the  higher  density  of  UN,  the  much  higher  thermal  conductivity  of  UN  over 
UOj,  and  the  expected  lower  fission  gas  swelling  with  bumup  of  UN.  The  higher 
density  of  UN  allows  for  higher  fission  densities  within  the  fuel  and  thus  a  smaller 
core  (and  mass)  for  producing  the  required  amount  of  power.  The  higher  thermal 
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Table  1 


SP-100  Specifications 


Number  of  Fuel  Elements 

3055 

Fuel 

UN 

Enrichment 

94% 

%  Theoretical  Density 

96% 

Cladding 

Nbl%Zr 

Liner 

Rhenium 

Fuel  Element  Length 

46.46  cm 

Fuel  Stack  Length 

7.62  cm 

Cladding  OD 

0.607  cm 

Cladding  Thickness 

0.055  cm 

Liner  Thickness 

0.010  cm 

Fuel  OD 

0.456  cm 

conductivity  of  UN  results  in  lower  temperature  gradients  in  the  fuel,  thus  allowing 
for  higher  operating  temperatures  (i.e.  higher  powers)  and  reducing  internal  thermal 
strain.  The  fission  gas  swelling  directly  affects  the  lifetime  of  the  reactor  and  UN’s 
lower  propensity  for  swelling  is  principally  derived  from  observations  from  the  early 
1970's.  Researchers  who  have  compared  UOj  with  UN  have  attributed  UN’s 
swelling  advantage  to  smaller  thermal  gradients  [2].  Higher  thermal  gradients  (and 
resulting  higher  fuel  temperatures)  in  UOj  introduce  aggravations  to  swelling  such 
as  lenticular  void  growth,  increased  fission  gas  and  bubble  mobilities,  and  phase 
changes.  Comparison  of  swelling  data  for  UC  with  that  of  UN  has  suggested  that,  at 
identical  temperatures,  UC  swells  at  a  rate  greater  than  both  UN  and  UOj  [2]. 

Fission  gas  swelling  is  potentially  the  largest  of  the  factors  which  contribute  to 
swelling  of  nuclear  fuel,  and  is  additionally  the  most  complex  and  difficult  to 
predict.  Unfortunately,  it  is  also  one  of  the  least  characterized  areas  of  UN. 

2.2  Boundary  Conditions 

The  first  task  in  performing  scoping  calculations  was  to  determine  the 
temperature  and  pressure  boundary  conditions  for  the  fuel  element.  The  computer 
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code  CENTAR  (£ode  for  Extended  Non-linear  Transient  Analysis  of 
Extraterrestrial  Reactor  System)  [3]  was  used  for  this  purpose.  CENTAR  was 
developed  to  model  space  reactor  systems  such  as  the  SP- 100.  However,  it  does  not 
contain  a  detailed  fuel  pin  model.  Instead  a  fuel  pin  is  treated  simply  as  a 
homogeneous  material  with  a  thermal  conductivity  and  pin  swelling  determined 
from  correlations,  which  are  based  on  existing  irradiation  data.  Because  of  this, 
their  accuracy  in  unusual  scenarios,  such  as  the  one  being  looked  at  here,  is 
unknown.  Fortunately,  the  boundary  conditions  as  calculated  in  a  systems  code  (i.e. 
coolant  temperature  and  pressure)  are  not  strong  functions  of  the  detailed  structure 
of  the  fuel  pins  and  thus  should  provide  reasonably  accurate  input  to  a  detailed  fuel 
modeling  code. 

Calculations  using  CENTAR  with  the  transient  scenario  considered  here 
indicated  that  at  10  kW  electric  (200  kW  thermal  and  the  minimum  needed  to  keep 
the  system  operational)  the  clad  surface  temperature  was  approximately  500°K, 
while  at  1  MW  (20  MW  thermal)  the  surface  temperature  jumped  to  approximately 
1500°K. 

2  J  Fuel  Temperature  Profile 

Using  the  calculated  outer  fuel  pin  temperatures,  the  steady-state  heat 
conduction  equation  was  solved  analytically  in  order  to  find  the  centerline 
temperatures  and  thermal  gradients  in  the  fuel.  The  thermal  conductivity  was  held 
constant  to  obtain  the  analytic  solution. 

The  temperature  increases  from  surface  to  centerline  (rf,^  -  were  5°K  at 
7*5=  500°K  and  300°K  at  7’^=  1500°K.  The  maximum  thermal  gradient  is  found  at 
the  fuel  surface.  These  values  are  38®K/cm  at  7’5=500°K  and  2550°K/cm  at 
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7’g=1500“K.  Based  on  these  numbers  one  would  expect  that  biased  effects  are  not 
significant  at  10  kW,  but  would  need  to  be  included  at  1  MW.  Since  the  low  power 
operation  will  determine  the  pre-transient  condition  of  the  fuel,  and  because  the 
transient  time  is  so  short,  the  following  discussion  will  concentrate  on  effects 
present  in  the  low  power,  low  temperature  phase  of  operation. 

2.4  Fission  Gas  Source 

In  order  to  model  a  transient,  the  pre-transient  condition  must  be  known. 

The  first  step  is  the  calculation  of  the  amount  of  fission  gas  deposited  in  the  grain  at 
200  kW  thermal  power  over  a  ten  year  life.  Dividing  200  kW  by  the  total  number  of 
fuel  pins  in  the  core  (3055),  an  average  thermal  rod  power  of  0.065  kW  results. 
Assuming  uniform  flux  within  the  fuel,  the  fission  rate  density  will  be  1.642x10*^ 
fissions/sec-cm^.  Using  a  grain  diameter  of  19  pm,  obtained  from  recent 
investigations  involving  SP-100  type  fuel  [4],  one  calculates  that  there  are 
185.7fissions/sec  within  the  grain.  Over  a  ten  year  period  this  translates  to 
5.850x  10‘°  fissions  in  a  grain.  Multiplying  this  by  a  stable  fission  gas  yield  of 
0.25fission  gas  atoms  per  fission  [5],  one  finds  that  1.462x10^°  fission  gas  atoms 
are  deposited  in  the  grain  over  ten  years,  representing  a  density  of  1.293x  10^  fission 
gas  atoms/cm^.  If  one  assumes  the  gas  atoms  to  be  evenly  dispersed  as  single  gas 
atoms  within  the  grain,  the  calculated  distance  between  two  gas  atoms  is 
19.8x  ICf’cm.  Comparing  this  to  the  diameter  of  a  xenon  atom  in  the  lattice 
(approximately  3x10'^  cm),  one  sees  that  local  effects  around  each  xenon  atom 
could  couple  with  other  xenon  atoms  which  could  be  highly  disruptive  to  the  lattice. 
If  disruption  is  not  expected  then  it  must  be  assumed  that  fission  gas  is  lost  from  the 
grain  or  that  bubbles  nucleate  and  grow  from  the  deposited  gas  atoms.  Since 
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disruption  of  the  lattice  is  not  expected,  diffusion  out  of  the  grain  and  nucleation 
must  be  considered. 

2.5  Steady-State  Diffusion  Out  of  the  Grain 

An  indication  of  how  much  gas  would  be  lost  from  the  grain  can  be  obtained 
by  applying  the  steady-state  diffusion  equation  to  the  single  atoms  deposited  in  the 
grain.  Taking  the  boundary  condition  for  diffusion  calculations  to  be  =  0  at  the 
grain  boundary,  one  finds  that  the  steady-state  number  of  fission  gas  atoms 
contained  within  the  grain  is  3.22xl(P^  atoms.  This  number  is  much  larger  than  the 
1.462x  10*°  fission  gas  atoms  deposited  over  ten  years,  so  it  is  clear  that  a 
steady-state  equilibrium  is  not  attained.  In  fact,  additional  random  walk  calculations 
show  that  each  xenon  atom  only  moves  a  net  distance  of  only  6.4x  10"^  pm,  far 
below  the  dimension  of  one  unit  cell  in  the  lattice.  Calculation  of  the  homogeneous 
nucleation  rate  based  on  the  fission  gas  concentration  and  thermal  diffusivity 
indicates  that  this  rate  does  not  significantly  alter  the  single  gas  atom  concentration. 
Ignoring  nucleation  for  the  moment  and  applying  the  time-dependent  form  of  the 
diffusion  equation  [6]  one  finds,  as  indicated  above,  that  there  is  no  significant 
release  of  fission  gas  to  the  grain  boundary.  Similarly,  calculation  of  the  thermal 
gradient  effects  show  that  biased  movement  of  fission  gas  atoms  is  insignificant  at 
this  temperature.  Scoping  calculations  thus  indicate  that  once  the  xenon  atoms  are 
deposited,  they  are  virtually  'frozen'  in  the  lattice  at  low  power.  However,  it  should 
be  noted  that  only  thermal  diffusion  of  fission  gas  was  considered  in  the  scoping 
calculations. 
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2.6  Heterogeneous  Deposition  and  Nucleation 

Although  manufactured  defects  can  act  as  heterogeneous  nucleation  sites,  the 
primary  source  of  defect  sites  will  be  the  deposited  atoms  and  fission  damage 
created  by  fission.  Accommodation  of  heterogeneous  effects  by  some  investigators 
has  been  done  by  including  semi-empirical  modifications  to  their  models  which 
account  for  fission  damage  and  heterogeneous  nucleation  [7].  Other  investigators 
have  justified  the  use  of  a  straight  homogeneous  approach.  Unlike  UO^,  UN  has  a 
greater  mobility  of  electrons  which  allow  the  energy  transferred  by  a  fission 
fragment  to  be  deposited  over  a  much  larger  region  [8,9].  As  a  result,  the  energy 
transferred  to  a  single  lattice  atom  is  lower  and  much  less  likely  to  impart  enough 
energy  to  dislodge  an  atom.  As  a  result,  lattice  damage  from  fission  will  be  lower, 
resulting  in  fewer  heterogeneous  nucleation  sites.  A  side  effect  of  this  property  is 
that  the  increased  electron  mobility  also  prevents  a  significant  thermal  spike  from 
being  produced  by  the  fission  fragment  As  a  result  the  resolution  of  bubbles 
should  occur  primarily  via  a  Nelson  knock-on  method  [  1 0]  versus  a  sputtering  or 
trapping  mechanism.  Because  of  this  difference  in  electrical  properties,  it  would  be 
questionable  to  simply  transfer  the  UOj  empirical  modifications  to  UN. 
Unfortunately,  modifications  tailored  to  UN  have  not  been  developed  due  to  the  lack 
of  appropriate  experimental  fission  damage  data  Because  fission  damage  is  lower 
in  UN,  heterogeneous  nucleation  will  occur  primarily  through  a  finite  probability 
that  a  deposited  gas  atom  comes  to  rest  next  to  an  existing  gas  atom  or  bubble. 
Calculations  of  this  probability  involve  the  ratio  of  nearest  neighbor  sites  for  a 
specific  bubble  size  to  the  total  number  of  lattice  sites  available.  Lattice  damage 
from  fission  is  not  expected  to  provide  a  significant  source  of  bubble  nucleation. 
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2.7  Transient  Regime 

Similar  calculations  to  those  already  discussed  were  performed  for  the 
transient  regime.  These  calculations  indicated  that  the  transient  time  at  high  power 
will  be  too  short  for  significant  movement  of  fission  gas.  Thus,  it  is  not  expected 
that  interactions  between  different  sized  bubbles  will  play  a  role  during  the  higher 
power  phase.  However,  the  predicted  amount  of  gas  release  is  dependent  on 
accurate  modeling  of  intragranular  fission  gas  behavior,  which  in  turn  depends  on 
several  properties  of  the  material.  Uncertainties  in  properties  need  to  be  considered 
in  fuel  design  to  determine  acceptable  safety  margins. 
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CHAPTER  III 

FISSION  GAS  MODELING  IN  UN 

3.1.  Historical  Overview 

Fission  gas  in  nuclear  fuels  has  been  of  interest  to  the  nuclear  community 
since  the  early  1960's.  The  theory  of  fission  gas  behavior  has  developed  from  a 
simple  model  of  single  gas  atoms  as  solid  fission  products  in  pseudo-amorphous  fuel 
[11,  12,  13]  to  the  current  state-of-the-art  of  detailed  models  describing  fission  gas 
behavior  both  within  fuel  grains  and  on  grain  edges  and  faces.  The  latter  models 
attempted  to  account  for  the  known  factors  which  affect  fission  gas  behavior,  i.e. 
resolution  from  bubbles,  heterogeneous  and  homogeneous  bubble  nucleation,  bubble 
coalescence,  bubble  interlinkage,  non-equilibrium  bubbles,  bubble  mobility, 
nonstoichiometric  effects,  grain  growth,  impurities  of  solid  fission  products,  bubble 
trapping  at  dislocations  and  defects,  fuel  cracking,  fission  enhancement  of  diffusion, 
etc.  Unfortunately,  the  complex  synergistic  nature  of  fission  gas  behavior  and  the 
difficulty  of  studying  the  fuel  at  a  microscopic  level  under  operating  conditions 
makes  identifying  and  describing  all  the  contributing  and  competing  variables  a 
monumental  task. 

3.1.1.  Modeling  Studies 

Although  some  of  the  existing  fuel  modeling  codes  are  based  on  empirical  or 
semi-empirical  correlations  to  describe  fission  gas  swelling  and  release,  these  are 
only  useful  within  the  range  of  operational  conditions  used  to  define  the  correlations 
and  lack  the  flexibility  to  perform  exploratory  analyses  in  new  areas  of  interest. 
Exploratory  analysis  requires  accurate  mechanistic  modeling  codes  based  as  much 
as  possible  on  known  theory.  It  is  for  this  reason  that  program  REDSTONE  was 
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designed  from  a  mechanistic  approach.  The  theoretical  basis  of  the  mechanistic 
equations  used  in  this  work,  and  for  the  majority  of  current  state-of-the-art 
mechanistic  fuel  modeling  codes,  dates  back  to  Gruber  [14]  who,  in  1966,  first 
assembled  a  set  of  equations  describing  diffusion  and  interactive  fission  gas 
processes  in  nuclear  fuel.  Subsequent  researchers  extended  and  adapted  Gruber’s 
technique  as  their  understanding  about  additional  fission  gas  processes  increased. 
Many  of  the  extensions  applied  by  subsequent  investigators  were  situation  specific 
and  thus  had  limited  usefulness  for  exploratory  analysis.  This  resulting  increase  in 
complexity  and  specificity  in  Gruber’s  method  appeared  to  continue  until  1979, 
when  emphasis  seemed  to  shift  from  extension  of  Gruber’s  method  to  simplification 
of  the  method  so  it  could  be  used  efficiently  in  large  integrated  fuel  modeling  codes. 
The  simplifications  were  primarily  of  two  types,  one  involving  collapsing  the 
individual  bubble  groups  into  a  single  average  bubble  group  and  the  other  involving 
utilizing  an  effective  diffusion  coefficient  for  single  gas  atom  calculations  to  account 
for  all  effects.  These  simplifications  were  developed  and  verified  using  calculations 
from  the  more  detailed  fuel  modeling  codes.  For  the  most  part,  all  were  based  on 
UOj  fuel  properties.  Both  Matthews  and  Wood  [15,16]  and  Rest  [17]  have  adapted 
Gruber’s  approach  to  their  own  codes,  TRAFIC  (TRAnsient  Fuel  Interpretative 
£ode)  and  GRASS  (^as  Release  And  ^welling  Subroutine).  Both  these  codes  have 
evolved  in  similar  directions  and  are  representative  of  the  mainstream  of  the  current 
mechanistic  approach  to  fission  gas  modeling.  The  processes  considered  in 
determining  the  effective  diffusion  coefficient  or  collapsing  of  bubble  groups  may  or 
may  not  be  the  same  in  non-U02  fuel  materials.  The  UO^  simplifications  might  only 
be  in  agreement  with  the  specific  cases  they  were  verified  against,  and  thus  would 
not  be  applicable  to  other  fuel  materials  and/or  operating  histories.  Before 
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simplifications  for  non-U02  fuel  can  be  applied,  initial  modeling  of  fission  gas 
behavior  is  best  done  by  reversion  to  the  more  detailed  equations  of  Gruber’s 
method  where  the  governing  equations  are  more  directly  related  to  experimentally 
measured  properties.  This  constraint  affects  the  direction  of  the  work  here  since  the 
fuel  of  choice  is  UN,  for  which  there  are  no  detailed  calculational  studies. 

3.1.2  Experimental  Studies 

As  with  the  modeling  developments,  experimental  studies  involving  nuclear 
fuels  have  concentrated  essentially  on  U02-  In  fact,  both  have,  out  of  necessity, 
gone  hand-in-hand.  As  a  result,  the  modeling  of  UO2  fuel  behavior  has  advanced  to 
a  point  that  is  quite  detailed,  but  which  cannot  be  duplicated  for  other  fuel  materials 
due  to  a  lack  of  appropriate  experimental  data.  Often,  material  properties  necessary 
for  applying  current  state-of-the-art  models  are  estimated,  usually  from 
corresponding  UOj  data  For  UN  there  is  some  support  for  substituting  UOj 
properties  as  some  researchers  have  suggested  that  certain  aspects  of  fission  gas 
behavior  in  UN  are  similar  to  those  of  UO2  [2,18].  In  addition  to  the  lack  of 
material  property  data  the  lack  of  experimental  studies  involving  non-UO.,  fuel 
materials  has  left  unresolved  what  processes  are  and  are  not  important  in  these 
materials.  Both  these  problems  have  to  be  dealt  with  in  the  case  of  UN.  The  small 
amount  of  UN  data  available  is  further  clouded  by  the  fact  that  early  investigators 
had  problems  in  producing  and  maintaining  pure,  stoichiometric  UN  [19,20,21], 

The  state  of  the  fuel  being  studied  was  not  always  well  characterized. 

3.1.3  Operational  Studies 

Studies  which  have  dealt  with  operating  scenarios  involving  power 
excursions,  in  addition  to  concentrating  on  UO2,  have  generally  considered  fuel 
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which  is  operating  initially  at  high  powers  (and  thus  high  temperatures)  and  then 
experiences  an  overpower  situation  to  very  high  temperatures  where  fuel 
restructuring  is  likely  to  occur.  The  controlling  processes  involved  in  these  studies 
are  different  than  those  in  the  investigation  here.  Micrographs  of  UN  fuel  irradiated 
at  relatively  low  temperatures  have  indicated  that  no  restructuring  of  UN  fuel  occurs 
at  low  temperatures  [19],  Predominantly  low  temperatures  are  expected  and  thus 
little  or  no  fuel  restructuring  is  expected  to  occur.  This  will  result  in  a  different 
pre-ramp  condition  of  the  fuel  than  in  other  studies.  Matthews  [15]  indicates  that 
the  form  of  the  intragranular  gas  at  the  start  of  a  transient  is  decisive  for  subsequent 
release  and  swelling.  No  studies  were  found  where  the  fuel  was  first  operated  for  an 
extended  period  at  the  low  temperatures  used  here  and  then  experienced  a  severe 
overpower  transient.  Even  if  such  studies  were  to  exist,  it  is  highly  probable  that 
they  would  have  concentrated  on  the  macroscopic  fission  gas  release  from  a  specific 
fuel  design,  probably  using  UOj,  and  not  on  the  microscopic  fission  gas  behavior  of 
UN  fuel. 

3.2  Fission  Gas  Behavior  Regimes 
3.2. 1  Intergranular  Fission  Gas 

For  intergranular  (grain  boundary)  fission  gas,  modeling  approaches  have  not 
congealed  into  any  single  mainstream  method,  possibly  due  to  the  greater  effect 
initial  manufacture  and  operating  conditions  have  on  open  porosity  and  grain 
boundary  dynamics.  Early  investigators  considered  the  intragranular  fission  gas  to 
be  released  once  it  reached  the  grain  boundary,  either  immediately  or  once  a 
specified  percentage  of  the  grain  boundary  had  been  filled  with  gas 
[22,23,24,25,26,27].  It  was  clear  from  accumulating  experimental  data,  however. 


18 


that  the  dynamics  of  fission  gas  behavior  on  grain  boundaries  were  much  more 
complicated  than  these  simple  approaches  suggested.  In  the  1970's,  Tucker 
[28,29,30]  looked  at  the  energies  involved  in  the  grain  boundary  dynamics  and  was 
able  to  outline  a  detailed  treatment  of  grain  face  and  grain  edge  bubbles.  Fission  gas 
release  occurred  only  after  grain  face  bubbles  had  interconnected  to  the  grain  edges 
and  grain  edge  bubbles  were  interconnected  to  the  exterior  or  central  void  of  the 
fuel.  His  approach  has  gained  some  acceptance  and  is  the  approach  taken  by 
Matthews  [31].  Unlike  the  dominance  of  Gruber’s  method,  however,  Tucker’s 
theories  are  not  universally  used  for  mechanistic  modeling.  At  least  one  investigator 
has  treated  the  grain  boundary  gas  using  difftisionai  rate  equations  similar  to  those 
used  in  Gruber’s  method  [32].  Since  intergranular  fission  gas  was  not  modeled  in 
the  study  here,  no  further  elaboration  of  intergranular  fission  gas  modeling  will  be 
made.  If  necessary,  the  amount  of  fission  gas  released  will  be  determined  using  a 
percentage  release  criteria 

3.2.2  Intragranular  Fission  Gas 

One  of  the  first  investigations  into  the  theory  of  inert  gas  behavior  in  nuclear 
fuels  was  performed  by  Greenwood  and  Speight  [33].  They  based  their  theory  on 
the  assumptions  that  bubbles  migrate  by  surface  diffusion,  there  was  no  resolution  of 
gas  from  bubbles,  the  gas  in  the  bubble  behaved  according  to  the  ideal  gas  law,  and 
the  bubbles  were  in  equilibrium  with  the  surrounding  matrix.  They  also  assumed  the 
material  to  be  a  perfect,  infinite  crystal,  and  that  the  bubble  distribution  could  be 
characterized  by  a  mean  bubble  radius.  Gruber  [14]  extended  this  work  to  eliminate 
the  depjendence  on  a  mean  bubble  radius  He  proposed  using  separate  rate  equations 
for  single  gas  atoms  and  each  bubble  size.  These  rate  equations  included  interaction 
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terms  between  single  gas  atoms  and  bubbles  to  produce  the  bubble  distribution. 

Other  assumptions  used  by  Greenwood  and  Speight  remained  assumptions  for 
Gruber’s  work.  Interaction  terms  included  by  Gruber  (nucleation  and  coalescence) 
were  based  on  colloidal  coagulation  work  performed  by  Chandrasekhar  [34].  The 
method  of  Gruber  provided  a  basis  that  described  the  detailed  bubble  distribution, 
yet  was  easily  adapted  to  specific  cases  by  simply  adding  additional  terms  to  the  rate 
equations.  The  specific  implementation  of  Gruber’s  method  used  here  will  be 
discussed  in  detail  in  the  chapter  covering  governing  processes. 

3^  Experimental  Data  on  UN 

3.3. 1  Overall  Status  of  UN  Property  Data 

Available  experimental  data  on  UN,  especially  for  the  specific  form  of  fuel 
proposed  for  the  SP-lOO,  is  limited  and  generally  not  sufficient  to  fully  characterize 
the  microscopic  processes  occurring  within  the  irradiated  material.  Advanced  fuels, 
such  as  UN,  were  first  considered  as  better  alternatives  to  UO,,  especially  for 
LMFBR  (Liquid  Metal  Fast  Breeder  Reactor)  technology  and  possible  applications 
in  space  [19,35,36,37],  in  the  late  1950’s  to  early  1960’s.  Early  investigations  of 
UN  and  UC,  as  well  as  the  mixed  fuels  (U,Pu)N,  and  (U,Pu)C  were  performed  in  the 
mid  to  late  1960's.  Unfortunately,  work  on  these  alternative  fuels  was  overwhelmed 
by  the  data  being  collected  on  UO,  and  (U,Pu)02  As  a  consequence,  investigations 
using  UN  dropped  off  in  the  early  1970’s  and  only  picked  up  again  in  the 
mid-1980’s  as  the  United  States  nuclear  space  program  and  SDI  (Strategic  Defense 
Initiative)  came  into  the  spotlight.  Most  of  the  experimental  data  on  UN  that  were 
obtained  in  the  1970’s  utilized  experimental  techniques  and  equipment  which  were 
not  as  refined  as  today,  and  the  effects  of  nonstoichiometry  and  fission  damage  in 
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UN  were  not  well  understood.  Unfortunately,  recent  investigations  utilizing  UN 
have  emphasized  macroscopic  observations  (i.e.  overall  release  and  pin  swelling)  of 
particular  fuel  designs  rather  than  the  microscopic  properties  which  would  be  useful 
for  generalized  fission  gas  modeling.  It  is  difficult,  if  not  impossible,  to  decipher  the 
microscopic  fission  gas  behavior  from  the  experimental  data  in  these  more  recent 
studies.  Overall,  experimental  data  available  on  UN  are  not  sufficient  for  building 
comprehensive  mechanistic  models  of  microscopic  fission  gas  behavior  in  UN  fuel. 
As  a  result,  some  of  the  data  necessary  must  be  borrowed  from  UO2  or  UC. 

Most  of  the  experimental  property  data  that  was  available  on  UN  up  to  1989 
has  been  compiled  into  a  database,  as  well  as  correlated  with  variables  such  as 
temperature  and  pressure,  at  Texas  A«&M  University  [38].  However,  some  of  the 
primary  data  needed  for  modeling  fission  gas  behavior  is  not  present  in  this 
compilation  and  had  to  be  gleaned  from  the  few  isolated  investigations  in  the 
literature.  Included  in  this  data  are  the  xenon  diffusion  coefficient,  the  surface 
diffusion  coefficient,  and  the  fission-enhanced  diffusion  coefficient.  No 
measurements  for  some  of  the  required  data  were  found.  Included  were  the  surface 
energy  in  UN,  the  fission  damage  dependence  of  properties  in  UN,  the  effect  of 
stoichiometry  on  UN  properties,  and  the  grain  boundary  energies  of  UN. 

Where  possible,  the  data  used  in  REDSTONE  was  taken  from  the  compiled 
database  at  Texas  Ai&M.  Other  data  were  taken  from  the  isolated  investigations  or 
taken  from  corresponding  UC  or  UOj  data. 

3.3.2  Diffusion  Coefficients 

The  diffusion  coefficients  of  U,  N,  and  the  fission  gases,  as  well  as  the 
surface  diffusion  coefficient  of  UN  are  absolutely  essential  to  modeling  fission  gas 
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behavior.  Each  of  the  different  types  of  diffusion  coefficients  are  discussed  in  the 
following  sub-sections. 

3.3.2. 1  Xe  Diffusion  Coefficient 

Data  for  fission  gas  diffusion,  whether  in  UN  or  other  fuels,  is  typically 
taken  to  be  the  diffusion  coefficient  of  xenon.  This  is  a  logical  choice  since  most  of 
the  stable  gaseous  fission  products  consist  of  Xe  (see  table  2).  Although  some  of  the 
short-lived  isotopes  have  half-lives  on  the  order  of  the  transient  times  considered 
here,  the  contribution  of  these  to  the  total  amount  of  gas  already  present  in  the  grain 
at  the  beginning  of  the  transient  will  be  negligible.  Unfortunately,  both  the  Xe 
diffusion  coefficient  and  the  surface  diffusion  coefficient  in  UN  are  poorly 
documented. 

Only  a  few  measurements  of  the  xenon  diffusion  coefficient  in  UN  have 
been  made.  Melehan  and  Gates  [20]  measured  the  in-pile  fission  gas  release  from 
sintered  uranium  nitride  in  1964.  Irradiations  and  gas  collection  were  performed  at 
several  temperatures  ranging  from  540°C  to  1700°C.  They  noted  that  the  fission  gas 
release  rate  up  to  800°C  was  independent  of  temperature.  This  was  attributed  to  a 
recoil  knockout  -  evaporation  release  mechanism.  Above  800”C,  the  release  rate 
followed  an  Arrhenius  behavior  described  by  the  equation 

£»xe  -4xl0~^exp|-'^  — j  cm^/s  (1) 

where  the  activation  energy  is  4.5xl0‘‘^  ergs.  Their  analysis  of  the  release  data  was 
based  on  the  theoretical  work  of  Beck  [39].  Although  an  equivalent  sphere  approach 
was  used,  the  size  of  the  spheres  was  determined  strictly  from  the  physical 
dimensions  of  the  fuel  element,  not  from  any  analysis  of  grain  size.  Also,  no 
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Table  2 


Isotopes  of  Krypton  and  Xenon  Released  in  Fission  [5] 


Isotope 

Half-Lifr 

Percent  Yield 

*5*Xe 

Stable 

3.2 

Stable 

4.7 

5.3  day 

6.6 

Stable 

6.6 

9.2  hr 

5.5 

‘“Xe 

Stable 

5.9 

Total  stable  xenon  yields 

20.4 

»^Kr 

Stable 

0.4 

*^Kr 

Stable 

0.85 

*^Kr 

Stable(10.6  yrs) 

0.15 

85m{^ 

4.4  hr 

1.3 

“Kr 

Stable 

1.4 

^Kr 

78  min 

2.5 

88Kr 

2.8  hr 

3.5 

Total  stable  krypton 

yields 

2.8 

correction  was  performed  to  account  for  interactions  of  gas  atoms  with  bubbles  or 
the  effect  from  nonstoichiometric  fuel.  Post-irradiation  examination  of  the  fuel  did 
not  reveal  evidence  of  density  changes  as  would  be  expected  from  decomposition  of 
UN.  However,  since  a  nitrogen  overpressure  was  not  maintained,  it  would  be 
expected  that  some  degradation  of  the  fuel  occurred  and  that  the  post-irradiation 
examination  was  not  able  to  detect  it.  Any  decomposition  would  affect  the  resulting 
calculated  diffusion  coefficient. 

About  the  same  time,  Biddle  [40]  measured  the  post-irradiation  emission  of 
xenon  from  uranium  carbonitride  powders  containing  varying  amounts  of  carbon. 
The  equivalent  sphere  model  developed  by  Booth  [41]  was  utilized,  with  the 
equivalent  sphere  size  being  determined  from  krypton  adsorption  measurements 
onto  the  surface  areas.  Biddle’s  analysis  did  not  account  for  gas  atom  interactions 
and  no  effort  appeared  to  be  made  to  keep  the  UN  stoichiometric  during 
measurement  Resulting  diffusion  equations  for  two  different  surface  areas  were 
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cm^  /  s;  Surface  area  =  1850  cm^  /  g  (  2  ) 
and 

£>xe  -  8  X 1 0~^  exp^  - ^  — j  cm^/s:  Surface  area  =  1000  cm^  /  g  .  (  3  ) 

These  two  results  appear  markedly  different.  However,  results  presented  for  all  the 
fuel  materials  seem  questionable  since  values  for  the  carbonitrides  exhibit  just  as 
wide  a  range  with  no  observable  trend  as  a  function  of  carbon  content. 

Measurements  of  xenon  diffusion  in  single  crystals  of  UN  were  performed 
by  Oi  [42].  The  evaluation  of  the  diffusion  coefficient  was  based  on  the  measured 
release  fraction  and  the  surface/volume  ratio  of  the  bulk  crystals.  No  allowance  was 
made  for  gas  atom  interaction  and  no  attempt  was  made  to  maintain  stoichiometric 
UN  during  measurements.  The  resulting  diffusion  coefficient  was 

£>xe  -  2. 05  X  10~^  exp|  -  — j  cm^  /  s  .  (4) 

This  equation  produces  values  higher  than  those  found  by  other  investigators.  This 
may  be  the  result  of  a  lack  of  traps  for  gas  atoms  in  a  single  crystal.  Oi  also 
mentions  that  the  scatter  in  their  data  was  considerable. 

Work  carried  out  in  1970  at  BMI  [43]  measured  fission  gas  release  from 
97  percent  dense  compacted  UN  which  had  been  neutron  activated  to  produce 
fissions.  The  resulting  data  were  analyzed  using  a  Hurst  [44]  model  (see  table  3). 
The  Hurst  model  approximates  the  effects  from  fission  gas  precipitation.  A  nitrogen 
overpressure  was  maintained  during  analysis  to  keep  the  UN  stoichiometric.  The 
diffusion  coefficient  for  UN  was  found  to  be 


£)xe  -  3  X 10  *  expl  - 


4.0x10 


-12  > 


kT 
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Table  3 

Estimated  Diffusion  Coefficients  for  UC,  UN,  and  UO., 


Fuel 

Estimated  Diffusion  Coefficient,  cm-/sec 

1400°C  1600“C 

UC 

3x10-'= 

8x10-" 

UN 

2x10-'^ 

2xl0-‘2 

uo. 

4x10-*^ 

2x10-" 

Z)xe  -  4.65 X  10  "‘exp 


4,9796x10 
kT 


-12\ 


cm^  /  s 


(5) 


Ritzman  [19]  notes  that  diffusion  coefficients  calculated  using  the  Hurst  model  are 
several  orders  of  magnitude  higher  than  values  previously  reported  due  to  the 
inclusion  of  precipitation  of  gas  atoms.  This  may  explain  the  lower  values  obtained 
by  Biddle  [40]. 

Later  diffusion  work  at  NASA  [45]  used  a  sweep  gas  facility  to  measure 
in-pile  krypton  release  rates  from  small  UN  specimens  irradiated  to  eight  percent 
bumup  and  at  high  fission  rate  density.  They  correlated  their  data  with  predictions 
of  a  fission  gas  release  model  which  included  atomic  diffusion,  fission-enhanced 
diffusion,  and  direct  recoil.  Calculation  of  any  one  of  these  contributors  depended 
on  the  accuracy  of  the  other  two.  Because  of  the  form  of  their  equation,  they  were 
required  to  calculate  the  effective  surface  area  for  release.  This  was  accomplished 
by  assuming  complete  interconnection  of  grain  boundary'  porosity.  No  correction 
was  made  for  gas  atom  interaction,  bubble  movement,  or  biased  motion,  even 
though  bubbles  were  clearly  evident  in  post-irradiation  analysis.  As  a  result,  the 
calculated  diffusion  equation  represents  overall  release  from  the  fuel  irregardless  of 
species.  The  diffusional  contributors  to  the  overall  gas  release  were  modeled  with 
the  steady-state  solution  to  the  diffusion  equation  in  rectangular  coordinates,  not 
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with  an  equivalent  sphere  model  as  used  previously.  The  resulting  diffusion 
equation  they  obtained  for  atomic  diffusion  was 


Dyj  -  5.6  X  10  ®exp  - 


1.08x10' 


cm^  /  s 


The  high  activation  energy  present  above  may  be  a  direct  result  of  neglecting 
contributions  from  bubbles  and  gas  atom  traps  such  as  grain  boundaries  (complete 
interconnection  was  assumed).  Also,  no  attempt  was  made  to  keep  the  UN 
stoichiometric,  and  the  small  size  of  the  specimens  may  have  caused  normally 
insignificant  processes  to  interfere  with  release.  In  addition,  it  should  be  noted  that 
the  diffusion  equation  presented  above  describes  krypton  diffusion,  not  that  of 
xenon.  Thus  some  difference  might  be  expected. 

Further  work  at  NASA  [46]  extended  the  above  analysis  to  include 
contributions  from  bubble  diffusion.  However,  other  criticisms  of  the  initial  study 
remained.  Their  adjusted  atomic  diffusion  equation  for  krypton  was 


Z)k,  -2.4xl0‘‘°exp  - 


2.60x10'^^ 


cm"  /  s 


It  is  unclear  if  the  inclusion  of  bubble  motion  was  enough  to  effect  the  large  changes 
observed  in  the  diffusion  equation.  The  latter  activation  energy  is  certainly 
consistent  with  those  of  other  investigators.  However,  the  pre-exponential  factor 
appears  to  be  far  too  low.  Differences  between  the  two  studies  were  not  discussed 
sufficiently  to  explain  the  change  in  diffusion  equations. 

The  dependence  of  diffusion  coefficients  on  the  amount  of  fission  damage 
present  has  been  used  as  one  explanation  of  the  differences  obtained  by  various 
investigators.  Experiments  with  UC  have  indicated  a  definite  dependence  of  the 
xenon  diffusion  coefficient  on  fission  damage  [47].  Unfortunately,  a  similar  trend 
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cannot  be  applied  to  UN  as  the  amount  of  fission  damage  in  the  preceding  study  was 
not  specified  and  fission  damage  work  with  UN  has  not  yet  been  performed.  Fission 
damage  dependence  is  not  included  in  program  REDSTONE. 

A  review  of  the  various  fission  gas  diffusion  equations  indicates  that  the 
latter  results  at  BMI  appear  to  be  the  most  reliable.  Not  only  were  steps  taken  to 
preserve  the  integrity  of  the  UN,  but  account  was  made  for  fission  gas  traps  not 
included  by  other  investigators.  Additionally,  the  form  of  the  fuel  used  at  BMI 
(sintered,  high  density,  stoichiometric  with  some  fission  damage  from  nc  it;  m 
activation)  was  similar  to  that  proposed  for  space  reactors.  The  diffusion  equation 
used  in  REDSTONE  is  that  derived  from  the  data  at  BMI. 

3. 3. 2.2  Irradiation  Enhanced  Diffusion 

Irradiation  enhanced  diffusion  of  fission  gas  in  nuclear  fuels  is  not  well 
characterized,  even  for  UOj.  Although  Dienes  and  Damask  [48]  have  looked  at  this 
problem  theoretically,  in  practice  irradiation  enhanced  diffusion  coefficients  are 
obtained  strictly  from  experiment.  In  these  cases,  irradiation  enhanced  diffusion  is 
used  to  explain  an  athermal  release  of  fission  gas  at  lower  temperatures.  For  UO.,, 
the  critical  temperature  for  irradiation  enhanced  diffusion  is  around  1000°C  [49,50], 
For  UN,  the  critical  temperature  appears  to  be  lower  around  800°C  [20].  Cornell 
[51]  derived  an  irradiation  enhanced  diffusion  equation  based  on  experimental 
observations  of  bubble  sizes  and  concentrations,  and  obtained  an  expression  that 
depended  on  the  fission  rate  and  resolution  rate.  The  irradiation  enhanced  diffusion 
coefficient  is  usually  expressed  in  terms  of  a  constant  and  the  fission  rate,  F. 

Typical  values  for  UO^  range  from  [49,52] 

De -2.0xl0-^^Fto2.0xl0-^‘F  cm^/s  (8) 
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where  Fis  the  fission  rate.  For  UN,  investigators  at  NASA  (National  Aeronautical 
and  Space  Administration)  [46]  appear  to  be  the  only  ones  to  have  calculated  an 
irradiation  enhanced  diffusion  rate.  The  value  they  found  for  UN  was 
Z>E -8.22xlO"^'F  cm^/s  (9) 

This  study  was  also  discussed  in  the  section  on  fission  gas  diffusion  which  involved 
measured  krypton  releases  using  sweep  gas  over  small  specimens  of  UN.  The 
results  were  fit  to  an  equation  which  included  contributions  from  irradiation 
enhanced  diffusion,  atomic  diffusion,  bubble  diffusion,  and  direct  recoil.  Although 
it  was  mentioned  that  these  contributions  were  interdependent,  evaluation  of  the 
irradiation  enhanced  diffusion  coefficient  was  accomplished  at  low  temperatures 
where  other  diffiisional  processes  were  negligible.  The  decrease  in  calculated 
irradiation  enhanced  diffusion  from  the  UO2  value  can  be  explained  from  the  lower 
amount  of  fission  damage  expected  in  UN  [49,53].  Since  no  other  source  for 
irradiation  enhanced  diffusion  in  UN  exists,  the  value  obtained  by  NASA  is  used  in 
REDSTONE. 

3.3. 2.3  Bubble  Diffusion  Coefficients 

Three  mechanisms  have  been  proposed  to  explain  bubble  mobilities:  vapor 
transport,  volume  diffusion,  and  surface  diffusion.  Vapor  transport  occurs  from  the 
evaporation  of  atoms  on  the  hot  side  of  a  bubble  and  deposition  on  the  cold  side. 
Since  surface  atoms  must  break  completely  away  from  the  surface,  high  activation 
energies  are  present,  and  thus  high  temperatures  are  required.  For  the  low 
temperatures  considered  here,  vapor  transport  should  be  negligible.  Experimental 
evidence  for  UOj,  where  the  driving  forces  and  temperatures  from  temperature 
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gradients  are  much  higher,  indicates  that  this  assumption  is  valid  up  to  moderately 
high  temperatures  [54,55], 

Volume  diffusion  occurs  from  the  bulk  movement  of  atoms  in  the  matrix 
around  the  bubble.  Calculations  on  the  relative  effects  between  surface  diffusion 
and  volume  diffusion  have  shown  that  only  when  bubbles  are  larger  than  about  1mm 
in  diameter  [15]  is  volume  diffusion  significant.  Since  only  small  bubbles  are 
expected,  volume  diffusion  is  not  considered  to  be  a  factor  in  the  study  here,  which 
leaves  only  surface  diffusion  to  explain  bubble  mobility.  By  far  the  most  widely 
accepted  model  for  bubble  mobility  is  the  surface  diffusion  mechanism.  Utilizing 
the  surface  diffusion  model,  the  bubble  diffusion  coefficient  is 
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(  10) 


where  =  bubble  diffusion  coefficient,  Z),  =  surface  diffusion  coefficient, 

A  =  atomic  spacing,  Q  =  atomic  volume,  and  =  bubble  radius. 

The  surface  diffusion  coefficient  for  UN  was  measured  by  Weaver  in  1968 
[56]  using  helium  implanted  UN.  After  annealing,  bubble  sizes  and  migration 
distances  up  a  thermal  gradient  were  measured  through  replication  electron 
microscopy  at  temperatures  ranging  from  1258°K  to  1858°K  and  thermal  gradients 
from  75°K/cm  to  880°K/cm.  Although  their  original  data  were  not  presented,  they 
did  calculate  a  temperature  dependent  surface  diffusion  coefficient  of 

2.93xl0-‘^) 


-  1.92  X  10  exp  -- 
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cm‘  /  s 
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Also  observed  was  that  coalescence  of  bubbles  did  not  occur  as  the  bubbles  traveled 
up  the  thermal  gradients.  Weaver  attributed  this  to  the  existence  of  large  stress 
fields  surrounding  nonequilibrium  bubbles.  This  latter  phenomena  has  been 
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theoretically  investigated  by  Ronchi  [57].  The  surface  diffusion  equation  obtained 
by  Weaver  is  used  in  program  REDSTONE. 

If  one  uses  the  above  equation  to  calculate  a  bubble  diffusion  coefficient  for 
a  bubble  the  size  of  one  xenon  atom,  one  gets  a  value  several  orders  of  magnitude 
higher  than  the  xenon  diffusion  coefficient.  Even  further,  Cornell  [51]  through 
observations  of  irradiated  fuel  postulated  that  very  small  bubbles  are  immobile  due 
to  being  trapped  at  dislocations,  and  thus  cannot  be  described  by  a  diffusion 
coefficient  as  given  above.  However,  Nichols  [58]  looked  at  much  of  the  data  on 
bubble  diffusion  and  concluded  that  small  bubbles  move  at  rates  far  greater  than 
those  assumed  by  Cornell,  though  still  not  at  rates  predicted  from  theory.  Lower 
small  bubble  mobility  is  probably  a  factor  since  small  bubbles  in  UO,  have  been 
observed  to  be  faceted  rather  than  spherical  [59,  60],  and  thus  would  not  be 
governed  by  the  same  enthalpies  as  those  of  larger,  purely  surface-diffusion 
controlled  bubbles.  This  presents  a  complication,  however,  since  it  is  general 
practice  to  treat  all  bubbles  as  being  controlled  by  the  same  mechanism.  So  that 
small  bubbles  are  accurately  modeled,  bubble  diffusion  coefficients  should  exhibit  a 
smooth  transition  from  single  xenon  atoms  to  the  larger,  surface  diffusion  controlled 
bubbles.  Various  theories  of  bubble  mobility,  some  based  on  other  diffusion 
mechanisms,  have  been  proposed  to  explain  experimental  observations  and  all  claim 
good  agreement  with  the  data  [61, 62].  No  predominant  theory  has  emerged. 

Gruber  [63]  and  Chien  [32]  have  adjusted  small  bubble  diffusion  coefficients  by 
considering  the  surface  atom  jump  distance  when  it  is  greater  than  the  circumference 
of  the  bubble.  Ronchi  [64]  adjusts  bubble  diffusion  coefficients  based  on  a 
transition  from  xenon  atom  to  a  bubble  size  large  enough  so  that  atomic  strains 
around  the  bubble  are  completely  relaxed.  Because  both  methods  give  similar 
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results,  Gruber’s  approach  is  applied  in  REDSTONE  due  to  its  easier 
implementation. 

Gruber’s  approach  was  to  apply  a  ‘rolling  stone’  model  for  surface  diffusion. 
Robertson  [65]  has  postulated  that  this  model  is  a  more  accurate  description  of 
surface  diffusion.  The  classic  model  of  surface  diffusion  is  that  it  occurs  one 
interatomic  distance  at  a  time  via  the  movement  of  vacancies  on  the  surface.  In  the 
‘rolling  stone’  model,  the  movement  of  surface  atoms  occurs  when  an  atom  pops  out 
of  the  surface  and  ‘rolls’  along  the  surface  to  an  existing  vacancy.  The  distance  of 
the  ‘roll’,  i.e.  the  mean  jump  distance,  will  span  10-1000  interatomic  distances. 
Gruber  shows  that  this  model  produces  results  consistent  with  existing  data. 

To  illustrate  the  approach,  a  bubble  of  radius  ‘/2„’  is  assumed  with  a  surface 
atom  mean  jump  distance  of  A,,  assumed  to  be  constant  and  independent  of  surface 
orientation  and  curvature.  The  net  movement  of  the  bubble,  A„,  can  be  calculated 
from  simple  geometry  from  the  arc  length.  A,,  and  the  bubble  radius,  /?„,  i.e. 

(12) 

For  the  above  becomes  A,,=A^.  For  A„  is  periodic.  Because  this 

result  yields  unrealistic  behavior,  Gruber  assumed  for  /?„<Aj/;rthat  A„=2/?„.  The 
correction  for  bubble  diffusion  can  then  be  represented  by 


d: 


2 


(13) 


31 


Substituting  the  appropriate  expressions  for  the  correction  factor  into  the  above 
equation  gives 


6AflA  A,  '1 

Jt 

(14) 
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The  only  additional  unknown  is  the  jump  distance,  Gruber  calculated  A,  by 
applying  his  modifications  to  the  data  of  Gulden  [61],  which  indicated  that  bubbles 
smaller  than  37  A  appeared  to  follow  an  dependence  rather  than  the 
dependence  expected  from  surface  diffusion.  Gruber  was  able  to  predict  the  trends 
in  Gulden’s  data  using  surface  diffusion  and  his  modifications.  He  was  also  able  to 
determine  A,  from  the  data  to  be  approximately  1 1 2  A .  Since  data  on  UN  to  perform 
a  similar  analysis  are  not  available,  it  is  assumed  in  REDSTONE  that  A^=  1 12  A  for 
UN. 

3. 3. 2. 4  Self-Diffusion  Coefficients 

Experimental  data  on  U  and  N  diffusion  in  UN  are  more  abundant  than  for 
fission  gas  diffusion.  The  available  data  have  been  surveyed  and  correlated  with 
nitrogen  partial  pressure,  temperature,  and  grain  size  in  a  database  recently  compiled 
at  Texas  A&M  University  [38].  Since  the  work  involved  in  compiling  the  database 
considered  data  quality  and  experiment  reliability,  it  was  decided  to  simply  use  the 
correlated  values  from  the  Texas  A&M  University  database.  An  extensive  review  of 
the  available  self-diffusion  data  will  not  be  presented  here.  Interested  readers  are 
referred  to  the  above  reference. 
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The  correlated  equations  for  self-diffusion  of  U  and  N  in  UN  are 
E\j  -  2.2146 — j  cm=/s  (16) 

and 

Z3^^  -  2.251  X  j  cm^/s  (17) 

where  P  =  nitrogen  partial  pressure  (atm),  7'=  temperature  (K),  and  g  =  average 
grain  size  (pm). 

3. 3. 2.5  Vacancy/Interstitial  Diffusion  Coefficients 

The  diffusion  of  vacancies  and  interstitials  is  poorly  documented  for  nuclear 
fuels.  Previous  investigators  adopted  values  from  other  materials.  Because  it  is  the 
vacancies  that  cause  bubble  growth,  several  investigators  who  have  considered  the 
effects  of  vacancies  have  not  included  interstitials.  Further,  only  a  few  investigators 
who  have  incorporated  vacancy  effects  have  published  the  actual  diffusion  data 
used.  Additionally,  vacancy  effects  are  often  incorporated  only  implicitly  by 
utilizing  creep  properties  of  the  fuel,  not  by  a  set  of  coupled  rate  equations.  As  a 
result,  there  is  slightly  more  documentation  concerning  vacancies  than  interstitials. 

For  interstitials,  both  Hayns  [66]  and  Dienes  and  Damask  [48]  use  the  same 
value  for  the  migration  energy,  3.20xl0’‘^  ergs,  where  Q'"is  the  energy  of 
migration  for  interstitials.  Hayns  took  his  value  from  austenitic  steel  while  Dienes 
and  Damask  took  theirs  from  a  representative  metal  ‘such  as  copper’  Since  the  two 
values  agree,  it  may  be  that  Dienes  and  Damask’s  representative  metal  was  not  all 
copper,  but  rather  a  conglomeration  of  data  available.  This  is  a  reasonable 
hypothesis  since  Dienes  and  Damask’s  work  predates  the  others  discussed  by 


33 


14years.  Hayns  used  a  pre-exponential  factor  of  10  while  Dienes  and  Damask 
used  a  theoretical  approach  to  arrive  at  10'^.  However,  several  values  needed  to 
arrive  at  the  latter  value  had  to  be  contrived.  Thus,  it  cannot  be  automatically 
assumed  that  the  theoretical  value  is  better  than  other  values  for  which  the  origin  is 
unknown.  Straalsund  and  Guthrie  [67]  assumed  a  pre-exponential  factor  of  one  for 
all  their  diffusion  equations  and  used  an  interstitial  migration  energy  value, 
reportedly  from  austenitic  steel,  which  was  different  from  the  above,  i.e. 

<^=2.Qn  xlO"*^  ergs.  The  reason  for  the  discrepancy  in  values  is  unknown. 
However,  they  mention  that  many  of  the  parameters  needed  are  highly  uncertain, 
and  thus  some  of  the  data  involves  reasonable  estimates.  Because  it  cannot  be 
determined  which  values  best  apply  to  UN,  and  given  the  agreement  between  Hayns, 
and  Dienes  and  Damask,  the  diffusion  of  interstitials  in  REDSTONE  was  chosen  to 
be  represented  by 


Dj  -  10  ^exp  - 
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None  of  the  above  works  which  considered  interstitials  were  concerned  with 
fission  gas  modeling.  Esteves,  et  al.  [68]  explicitly  considered  vacancies  (but  not 
interstitials)  in  modeling  fission  gas  behavior  in  UO,.  Their  form  of  the  diffusion 
equation  for  vacancies  was  derived  by  taking  the  pre-exponential  factor  from 
uranium  self-diffusion  and  a  vacancy  migration  energy  equal  to  the  vacancy 
formation  energy,  i.e. 


ZJy  -  4  X  10"’exp(  - 
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where  0^=0^  =  2.44xl0'‘^  ergs.  The  energy  of  formation  was  taken  as  one-half  the 
activation  energy  for  self-diffusion.  The  basis  for  their  pre-exponential  factor  is 
questionable  since  Olander  [5]  notes  that 

Z),  -  D^C^Q  .  (  20  ) 

Thus,  one  would  not  expect  the  self-diffusion  pre-exponential  factor  and  the  vacancy 
pre-exponential  factor  to  be  the  same.  Additionally,  it  would  be  expected,  as  is  the 
case  with  the  other  investigators  reviewed,  that  the  migration  energy  would  be 
approximately  80%  to  90%  of  the  formation  energy.  Other  investigators  took  their 
values  from  other  metals.  Dienes  and  Damask  [48]  again  used  values  from  a 
‘representative  metal’.  Their  theoretically  based  pre-exponential  factor  was  one  and 
the  vacancy  migration  was  1.28x10'*^  ergs.  As  will  be  seen,  this  value  is  much 
lower  than  others.  Because  of  the  previously  mentioned  uncertainties  in  the  source 
of  Dienes  and  Damask’s  data,  their  values  for  vacancy  diffusion  were  not  considered 
for  use  in  REDSTONE.  Straalsund  and  Guthrie  [67],  Hayns  [66],  and  Brailsford 
and  Bullough  [69]  all  use  values  derived  from  steel.  Straalsund  and  Guthrie 
assumed  a  pre-exponential  factor  of  one,  and  used  a  migration  energy  of 
Q!^2.Q5  xlO'*^  ergs.  Hayns,  and  Brailsford  and  Bullough  both  used  a 
pre-exponential  factor  of  0.6.  However  their  migration  energies  were  slightly 
different  with  Hayns  using  2.24x10"'^  ergs,  and  Brailsford  and  Bullough  using 
2^?*  =  2.08x10"*^  ergs.  The  discrepancy  between  these  last  three  migration  energies 
could  not  be  resolved,  although  it  is  possible  that  the  value  used  by  Straalsund  and 
Guthrie  involves  a  ‘reasonable  estimate’. 
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The  vacancy  migration  energy  used  in  REDSTONE  is  an  average  of  these 
three  energies.  The  pre-exponential  factor  was  taken  as  0.6,  i.e. 

By  — j  cm^/s.  (21) 


3.3.3  Miscellaneous  Data 

Non-diffusional  experimental  data  is  also  required  to  adequately  model 
fission  gas  behavior.  The  more  important  ones  are  discussed  below. 


3.3.3. 1  Thermal  Equilibrium  Concentrations 

The  equilibrium  concentrations  of  vacancies  and  interstitials,  although 
insignificant  in  an  irradiation  environment,  are  needed  to  specify  the  boundary 
conditions  of  bubbles,  dislocations,  and  grain  boundaries.  The  equilibrium 
concentrations  are  governed  through  thermal  emission  via 

and 


CT-;^exp( 


(23) 


where  Q  =  atomic  volume,  =  vacancy  energy  of  formation,  and  Qj  =  interstitial 

energy  of  formation.  As  with  diffusion  data,  the  energies  of  formation  are  poorly 
documented.  This  is  particularly  true  for  interstitials.  Olander  [5]  notes  that,  unlike 
vacancy  formation,  a  variety  of  interstitial  sites  are  available  and  thus  it  is  harder  to 
determine  the  interstitial  energy  of  formation.  Again,  investigators  have  used  data 
available  from  similar  metals.  Straalsund  and  Guthrie  [67]  used  a  value  from 
austenitic  steel,  6.9x10"'^  ergs,  while  Dienes  and  Damask  [48]  used  a  value 
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from  their ‘representative  metal’,  6.4xl0''2  ergs.  Olander  notes  further  that  the 
energy  of  formation  for  interstitials  is  probably  higher  than  that  for  vacancies  in  all 
solids.  In  practice  for  nuclear  fuels,  the  equilibrium  concentration  of  interstitials  is 
taken  as  zero  due  to  the  high  formation  energy  [5,48,70,71].  This  same  practice  is 
used  in  REDSTONE,  i.e.  Cf  =  0. 

For  vacancies,  the  equilibrium  concentration  is  retained.  Hayns  [66]  and 
Brailsford  and  Bullough  [69]  used  the  same  values  from  steel,  Ql  =  2.56x10'*^  ergs. 
Straalsund  and  Guthrie  [67]  used  energies  of  formation  in  the  range  Ql  =  2.6x1 0"*^ 
to  2.8x10"*^  ergs.  Dienes  and  Damask  [48]  used  values  from  their  ‘representative 
metal’,  ^  =  1.92x10"*^  ergs.  Like  their  value  of  the  vacancy  migration  energy,  this 
value  is  much  lower  than  that  used  by  other  investigators.  Because  of  this 
discrepancy,  this  value  was  not  considered  for  use  in  REDSTONE.  As  discussed 
earlier,  Esteves,  et  al.  [68]  took  half  of  the  activation  energy  for  self-diffusion  in 
UOj  as  the  formation  energy  for  vacancies  [72],  i.e.  ^  =  2.44xl0'*2  ergs.  Another 
fission  gas  modeling  work.  Ting,  et  al.  [73],  modeled  UC  using  a  vacancy  formation 
energy  provided  by  DonnerandSchiile[  74],Qj  =  2.71  x lO  '-  ergs. 

All  the  values  above,  except  for  that  of  Dienes  and  Damask,  fall  in  the  same 
general  range.  Since  one  value  does  not  appear  any  more  applicable  than  another, 
the  vacancy  energy  of  formation  used  in  REDSTONE  is  an  average  of  the  above 
values,  with  the  exception  of  Dienes  and  Damask’s  value,  i.e.  =  2.60xl0'‘^  ergs, 
with 
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3.3. 3.2  Surface  Energy 

The  surface  energy  is  important  for  determining  bubble  sizes  and  excess 
pressures.  No  available  data  on  the  surface  energy  of  UN  were  found.  Surface 
energy  data  is  available  concerning  the  other  two  nuclear  fuels,  UO,  and  UC.  An 
early  experimental  study  [75]  found  the  free  surface  energy  of  UOj  to  be 
640ergs/cm  ^  and  the  free  surface  energy  of  UC  to  be  1000  ergs/cm^.  These 
numbers  are  very  approximate,  with  stated  uncertainties  on  the  order  of  20  to  30 
percent.  The  free  surface  energy  of  UO2  was  calculated  from  theory  by  Benson  [76] 
to  be  1030  erg/cm^.  Warner  and  Nichols  [77]  estimated  the  surface  energy  to  be 
approximately  1000  ergs/cm^.  About  the  same  time,  Eberhart  [78]  found  critical 
surface  energies  of  UOj  to  lie  between  600  and  1650  ergs/cm^  depending  on 
whether  an  argon  or  hydrogen  atmosphere  was  present.  Later,  Maiya  [79] 
performed  experiments  of  the  growth  of  grain  boundary  grooves  and  determined  an 
average  free  surface  energy  of  626  ergs/cm^.  Clearly,  a  consistent  agreement  on  the 
value  of  free  surface  energy  in  UOj  has  not  been  reached.  Most  often,  the  value 
used  in  fuel  modeling  codes  is  1000  ergs/cm“.  For  UC,  reported  surface  energy 
values  have  varied  just  as  much.  Initial  studies  set  the  free  surface  energy  at 
730ergs/cm  ^  [80]  to  780  ergs/cm^  [81].  However,  further  work  indicated  that  the 
free  surface  energy  of  UC  was  as  high  as  2000  erg/cm‘  [82,83].  Although  some  fuel 
modeling  codes  have  used  the  lower  values  for  UC,  the  spread  of  available  data 
suggests  that  the  use  of  1000  ergs/cm^  would  not  be  unreasonable.  Since  no 
reported  data  for  free  surface  energy  of  UN  have  been  found,  and  consistent  values 
are  not  available  for  either  UOj  or  UC,  the  value  for  free  surface  energy  of  UN  used 
in  REDSTONE  is  1000  ergs/cm^. 
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3. 3. 3.3  Lattice  Parameter 

As  with  self-diffusion,  the  lattice  parameter  of  UN  has  also  been  correlated 
with  temperature  at  Texas  A&M  University  [38].  The  equation 
a-4.879  +  3.264xl0”’7’  +  6.889xl0-V^  A  (25) 

is  used  in  REDSTONE  where  T  =  temperature  (K). 

3. 3. 3. 4  Yield  Stress 

The  yield  stress  is  used  in  REDSTONE  to  determine  the  minimum  sizes  of 
bubbles.  The  yield  stress  for  UN  was  approximated  from  the  stress-strain  curves  of 
Werner  and  Blank  [84].  Values  consistent  with  the  plastic  behavior  were  read  from 
the  curves  and  then  correlated  with  a  least  squares  fit.  The  resulting  correlation  is 
ay --9.256x10*7  +  1.735x10*2  Pa.  (26) 

The  values  produced  with  such  a  correlation  are  very  crude.  Fortunately,  bubble 
sizes  change  very  little  with  large  changes  in  pressure  near  this  stress.  Thus, 
inaccuracies  in  the  yield  stress  should  have  minimal  effect  on  the  program  results. 
The  main  effect  of  the  yield  stress  is  to  keep  the  bubble  sizes  from  becoming  too 
small. 

3. 3. 3. 5  Thermal  Conductivity 

The  thermal  conductivity  of  UN  has  been  correlated  with  temperature  and 
porosity  at  Texas  A&M  University  [38].  This  correlation  was  used  in  determining 
the  temperature  profile  in  the  fuel  during  scoping  calculations.  The  thermal 
conductivity  of  porous  UN  is  given  by 

itp -1.864r®“®‘exp(-2.147)  W/m-k  (27) 

where  7  =  temperature  (K)  and  P  =  volume  fraction  porosity. 
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CHAPTER  IV 

GOVERNING  EQUATIONS  AND  TERMS 


4.1  Governing  Processes 

The  equation  which  is  the  basis  of  program  REDSTONE,  the  diffusion 
equation,  can  be  considered  a  balance  of  production  and  loss  rates  of  individual 
species,  be  they  gas  atoms,  gas  bubbles,  vacancies,  or  interstitials.  Considering  first 
individual  gas  atoms,  a  balance  of  production  and  losses  within  the  grain  typically 
would  consider  the  following; 
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For  UN  fuel  and  the  operating  scenario  to  be  considered  here,  some  of  the 
above  processes  can  be  eliminated  from  consideration.  Based  on  scoping 
calculations,  the  amount  of  intergranular  gas  would  be  so  small  that  resolution 
effects  from  these  bubbles  would  be  negligible.  Assuming  no  grain  boundary 
resolution  also  provides  for  the  worst  case  situation  where  diffusion  out  of  the  grain 
is  maximized.  Within  the  UN  fuel  gram,  the  effect  of  crystalline  defects  on  gas 
atoms  will  be  small,  in  the  case  here  assumed  negligible,  due  to  the  lower  amount  of 
fission  damage  expected  in  UN  [9].  The  dominance  of  fission  gas  bubbles  will 
provide  for  the  major  capture  defect  [85].  The  effect  of  fission  damage  defects  is 
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somewhat  accounted  for  as  the  Xe  diffusion  coefficient  used  was  measured  in 
irradiated,  fission  damaged  material  [43],  Loss  by  grain  growth  can  also  be 
eliminated,  since  the  temperatures  of  interest  during  the  long  term  operation  are  too 
low  for  significant  grain  growth  to  occur  [20,37],  During  the  short  term  at  high 
power,  the  time  is  too  short  for  any  significant  grain  growth.  After  eliminating  the 
above  terms,  a  simplified  equation  can  be  written: 
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(29) 


Biased  effects  are  implicitly  included  in  the  above.  Modeling  of  biased  migration 
effects  can  likewise  be  simplified  since  biased  migration  would  be  minimal  in  UN  at 
low  temperatures  and  at  high  temperatures  the  time  available  severely  limits  atom 
and  bubble  movement.  Specific  simplifications  to  biased  effects  will  be  discussed  in 
sections  dealing  with  those  effects. 

For  gas  bubbles,  a  similar  production- loss  equation  can  be  written: 
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The  same  arguments,  as  for  single  gas  atoms,  can  be  used  to  eliminate 


effects  from  crystalline  defects  and  grain  growth.  The  bubble  rate  equation  can  thus 


be  simplified: 
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(31) 


Vacancies  and  interstitials  can  be  treated  likewise  with  an  overall  equation: 
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Unlike  gas  atoms  and  bubbles,  effects  from  crystalline  defects  are  included  as  they 
provide  the  driving  force  for  bubble  growth  [5,86].  Matthews  and  Wood  [15]  note 
that  the  effect  of  crystalline  defects  as  vacancy  sources  will  be  negligible  compared 
to  the  grain  boundary  for  grains  less  than  lOOpm.  However,  vacancy  effects  from 
crystalline  effects  are  important  for  determining  the  vacancy/interstitial  balance,  and 
thus  vacancy  emission  from  defects  is  included.  Release  of  interstitials  from  sinks 
are  considered  negligible  since  the  equilibrium  concentration  of  interstitials  in  UO,, 
and  presumably  UN,  is  negligible  [5,87].  Release  of  vacancies  from  bubbles  and 
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grain  boundaries  are  not  incorporated  explicitly,  but  as  boundary  conditions  when 
solving  equations  for  the  overall  concentrations.  Thus,  the  simplified  equation  for 
vacancies  and  interstitials  is; 
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The  rate  equation  describing  bubble  vacancy  adjustment  is  a  modified 
version  of  the  simplified  bubble  rate  equation; 
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(34) 


The  actual  number  of  vacancies  in  a  single  bubble  is  obtained  by  dividing  the 
solution  to  the  above  equation  by  the  bubble  concentration. 

The  form  of  the  equation  used  in  REDSTONE  is  the  calculation  of  bubble 
radii,  not  the  calculation  of  bubble  vacancies.  The  bubble  radius  depends  directly  on 
the  number  of  vacancies.  The  radius  equation  can  thus  be  related  to  the  above 
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equation  describing  bubble  vacancies.  Many  of  the  properties  discussed  in  earlier 
chapters  also  depend  on  the  bubble  radius.  The  bubble  radius,  in  turn,  indirectly 
depends  on  the  fuel  properties.  Calculating  the  bubble  radius  directly  allows  one  to 
update  these  properties  concurrently,  improving  the  stability  of  calculation.  The 
specifics  of  converting  from  the  bubble  vacancy  equation  to  the  radius  equation  will 
be  covered  in  later  sections  and  chapters. 

4.2  1-D  Diffusion  Equation  in  Spherical  Coordinates 

The  diffusion  equation  for  a  single  species  (n=l  for  atoms,  n>l  for  bubbles, 
n=v,i  for  vacancies,  interstitials  respectively)  can  be  written  in  general  terms  as 

N  N  N  N 

-  v-(  D.vcJ  +  a.  +  I  +  II 

Ica^Q  ki«n  ka*nj#n 

(35) 

where  C„  =  the  unknown  quantity  of  n’th  species  (concentration  or  vacancies  per 
bubble),  £>„  =  diffusion  coefficient,  =  constant  factors,  /3„  =  constant  multipliers 
of  the  unknown  quantity,  =  constant  multipliers  of  the  square  of  the  unknown 
quantity,  =  constant  multipliers  of  products  of  the  unknown  quantity  and  k’th 

species,  =  constant  multipliers  of  products  of  k’th  and  j’th  species,  =  constant 

multipliers  involving  the  k’th  species.  Note  that  the  ‘j’  index  in  the  double 
summation  term  begins  at  ‘k’  rather  than  ‘1’.  Since  there  is  no  distinction  between 
contributions  from  ‘C^Cj’  or  ‘CjC^’,  starting  the  summation  at  ‘  1’  would  wrongly 
include  a  second  contribution  from  this  term. 
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For  1-D  spherical  coordinates  the  gradient  term,  which  describes  the 
diffusion  loss  out  of  the  grain,  becomes 


(36) 


Since  the  diffusion  coefficient  is  held  constant  during  calculations,  this  equation 
becomes 


V  •  (Z)  VC  )  -  ^ p-  — 

^  P  dr\  dr  I 


(37) 


The  specific  multipliers  of  the  overall  equation  are  determined  by  the 
individual  processes.  These  processes  are  listed  in  table  4  by  species  type  and  which 
multiplier  is  affected.  Individual  processes  are  discussed  in  more  detail  by  equation 
type  below. 


43  Gas  Atom  and  Bubble  Concentrations 

With  the  exception  of  biased  diffusion  out  of  the  grain,  the  processes 
governing  the  change  in  gas  atom  and  bubble  concentrations  are  described  by 


^  -  D„V"C„  -  K,C„  H. 

Ol 

k«n 


(38) 


k-n 


where  each  of  the  terms  are  in  order  of  appearance:  diffusion,  source  loss  from 
producing  the  next  larger  bubble,  source  gain  from  the  next  smaller  bubble, 
resolution  loss  to  the  next  smaller  bubble,  resolution  gain  from  the  next  larger 
bubble,  self-coalescence  loss,  coalescence  loss  from  random  and  biased  effects,  and 
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Table  4 


Identification  of  Processes 
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Bubble  V  ic  ;cy 
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a 
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BUBBLES 

coalescence  gain  from  random  and  biased  effects.  Each  will  be  discussed  in  the 
following  sections. 

4.3. 1  Production  by  Fission 

The  production  of  gas  atoms  by  fission  can  be  represented  by 
K  -  fpo^  (  39  ) 

where  ^FG  =  yield  of  fission  gas  atoms  per  fission  and  F-  fission  rate.  Both  the 
yield  and  fission  rate  are  taken  as  constant  within  a  time  step  and  spatially.  This 
latter  constraint  is  consistent  with  the  fast  neutron  flux  expected  in  the  reactor  core. 
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The  yield  of  fission  gas  atoms  per  fission  was  taken  as  0.25  [5]  from  consideration 
of  the  stable  isotopes  of  kiypton  and  xenon. 

Researchers  generally  have  considered  the  fission  gas  source  to  contribute 
only  to  single  gas  atoms.  However,  scoping  calculations  revealed  that  if  gas  atoms 
are  immobile,  then  the  final  density  of  atoms  is  too  great  to  be  realistically  supported 
by  the  lattice.  Because  of  the  high  density,  there  is  a  distinct  possibility  that 
deposition  of  gas  atoms  in  latter  stages  of  operation  would  have  a  finite  probability 
of  forming  new  bubbles  or  increasing  the  size  of  existing  bubbles.  Such  a  process 
can  be  modeled  by  considering  the  overall  fission  gas  source  to  be  divided  up  into 
separate  components  for  single  atoms  and  each  bubble  size.  Apportionment  is 
determined  by  calculating  the  ratio  of  nearest  neighbor  lattice  sites  for  a  particular 
bubble  size  to  the  total  number  of  lattice  sites  available  within  a  unit  volume. 

Lattice  sites,  whether  filled  or  unfilled,  are  assumed  to  always  exist  and  no 
calculation  is  done  to  alter  the  dimensions  of  the  fuel  grain.  Since  the  number  of 
lattice  sites  remains  constant,  the  total  number  of  lattice  sites  in  a  unit  volume  can  be 
calcultited  once  at  the  beginning  of  calculations,  i.e. 


where  is  the  lattice  constant  The  number  of  nearest  neighbor  lattice  sites  for  a 
single  gas  atom  in  a  face-centered-cubic  (FCC)  structure  is  12  [5].  The  number  of 
nearest  neighbor  lattice  sites  for  a  bubble  can  be  calculated  from  the  surface  area  of 
the  bubble  and  the  (average)  density  of  lattice  sites  on  the  surface.  When  multiplied 
by  the  bubble  concentration,  C„,  this  gives  the  total  nearest  neighbor  sites  for  a 
bubble  size  in  a  unit  volume.  The  bubble  surface  area  is 


\ 


(41) 
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where  is  the  radius  of  bubble  with  ‘n’  atoms.  In  calculating  the  density  of  lattice 
sites  on  a  bubble  surface,  most  investigators  consider  only  the  (100)  crystal 
direction.  A  more  accurate  density  can  be  calculated  by  considering  the  three 
principal  directions  (100),  (110),  and  (111)  and  taking  an  average.  Lattice  site 
surface  densities  for  each  direction  are; 


2_  2 


(42) 


(100)  (110)  (111) 

The  average  of  these  three  directions  is  2. 10032  a^.  Combining  with  the  bubble 
surface  area,  the  number  of  nearest  neighbor  sites  to  a  single  bubble  is  given  by 
Nearest  Neighbors  -  8.401 28jr^ao^  .  (  43  ) 

Multiplying  by  the  bubble  concentration  and  fission  gas  source,  and  dividing  by  the 
total  number  of  lattice  sites,  the  source  term  for  bubble  one  atom  larger  than  ‘n’  is 
given  by 


-finCn  -  8.40128;i/?Ja,C„rpoF 


(44) 


Correspondingly,  the  source  for  single  atoms  deposited  in  the  grain  is  given  by 
A:„  -  a,  -  -  YraF- 28.40i:8jifl;o.C„rpo^ 


n«l 


n«l 


(  45) 


4.3.2  Nucleation 

Historically,  nucleation  of  bubbles  has  been  treated  homogeneously  with 
modifying  terms  added  to  account  for  heterogeneous  effects  [88,  89,  90]. 
Heterogeneous  nucleation  has  periodically  been  explicitly  applied  in  an  effort  to 
explain  anomalies  in  observed  trends  [91].  Heterogeneous  nucleation  was  utilized  in 
the  BUBL  code  [87].  However,  this  code  also  assumed  gas  release  from  the  grain 
occurred  entirely  from  bubble  motion  and  bubble  interactions  occurred  only  at 
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dislocations.  Single  gas  atoms  in  the  matrix  only  served  to  nucleate  and  increase  the 
size  of  bubbles  on  dislocation  lines.  Robertson  [92]  makes  the  observation  that 
applied  heterogeneous  nucleation  is  on  a  scale  so  fine  that  little  difference  exists 
between  the  two  approaches  in  practice.  This  is  supported  by  the  calculations  of 
Dollins  [93]  which  indicated  that  a  gas  atom  is  several  orders  of  magnitude  more 
likely  to  meet  another  gas  atom  than  a  dislocation.  Also,  Warner  and  Nichols  [94] 
note  that  the  parameter- fitting  used  to  adjust  models  to  actual  experimental  data 
minimizes  any  large  discrepancies  between  the  two  models.  In  REDSTONE, 
homogeneous  nucleation  is  applied  with  the  only  heterogeneous  effect  being  the 
probabilistic  production  discussed  in  the  previous  section.  This  assumption  appears 
legitimate  since  fission  damage  will  be  lower  in  UN  than  in  UOj  and  crystalline 
defect  concentration  will  be  negligible  compared  to  the  bubble  density. 

Typically,  a  stable  nucleus  of  a  bubble  is  considered  to  be  a  diatomic  cluster 
of  gas  atoms  [88].  All  researched  modeling  efforts  have  assumed  this,  and  this  will 
be  assumed  to  be  the  case  here.  The  modeling  of  bubble  nucleation  will  introduce  a 
loss  term  into  the  single  gas  atom  equation  and  a  gain  term  into  the  equation 
describing  diatomic  gas  bubbles.  Gain  of  diatomic  bubbles  through  nucleation  is 
described  by  an  equation  given  by 

Gain  Nucleation  Rate  =  k\iC^  (  46  ) 

where  C,  =  gas  concentration  and  ,  =  rate  constant.  The  nucleation  loss  term  for 
single  gas  atoms  is 

Loss  Nucleation  Rate  -  -2it,  .  (  47  ) 

The  equation  is  multiplied  by  two  because  two  gas  atoms  are  removed  for  every 
diatomic  cluster  formed. 
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Following  a  purely  theoretical  approach,  Olander  [5]  derived  the  rate 
constant  for  single  vacancy  processes,  in  terms  of  a  combinatorial  number  and 

basic  properties  of  the  fuel  and  gas  atom,  i.e. 


where  =  diffusion  coefficient  of  vacancies,  z^  =  a.  to-be-determined 
combinatorial  number,  Q  =  volume  per  lattice  site  =  a^/4  for  a  FCC  lattice  structure, 
and  =  lattice  parameter.  The  combinatorial  number  as  calculated  by  Olander  is 
based  on  the  number  of  nearest  neighbors,  the  species  jump  frequency,  and  the 
species  concentration.  For  a  FCC  structure,  the  calculated  combinatorial  number  is 
84.  In  considering  interactions  between  xenon  atoms,  Olander  states  that  the 
combinatorial  number  is  probably  higher  than  an  atomic  value  (such  as  for 
vacancies)  because  evidence  suggests  that  the  xenon  atoms  migrate  as  a  complex 
involving  the  xenon  atom  and  several  vacancies.  Researchers  have  estimated  that  a 
single  gas  atom  is  combined  with  two  to  four  vacancies,  with  the  most  probable 
being  a  neutral  trivacancy  in  UOj  [95,96].  For  ON,  there  have  not  been  any  studies 
to  determine  the  optimum  number  of  vacancies  associated  with  a  gas  atom.  Since  a 
neutral  set  of  vacancies  seems  most  logical,  a  gas  atom  in  UN  will  be  assumed  to 
consist  of  the  atom  and  a  neutral  divacancy.  Olander  does  not  calculate  the 
corresponding  combinatorial  number  for  a  xenon-vacancy  complex  nor  does  he  give 
an  indication  of  how  much  larger  this  number  should  be  than  the  atomic  value.  If 
the  number  of  occupied  nearest  neighbors  to  the  complex  are  the  controlling  factor, 
as  indicated  by  Olander’s  analysis,  then  a  divacancy  complex  in  UN  would  have  18 
versus  12  nearest  neighbor  sites.  Based  on  this  change  alone,  the  combinatorial 
number  becomes  126.  However,  this  is  an  overly  simple  adjustment  and  there  are 
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probably  other  mitigating  factors,  such  as  the  effect  of  divacancies  associated  with 
the  jumping  atom  as  well  as  whether  a  gas  atom  jumping  into  a  nearest  neighbor  to  a 
vacancy  constitutes  nucleation,  which  might  reduce  the  combinatorial  number. 

Other  researchers  have  successfully  utilized  a  single  vacancy  mechanism 
(combinatorial  number  =  84)  to  describe  xenon  migration  [88].  Because  of  the  lack 
of  evidence  concerning  xenon-vacancy  complexes  in  UN  and  thus  the  calculation  of 
a  combinatorial  number,  a  combinatorial  of  84  will  be  assumed.  In  order  to  convert 
from  the  vacancy  form  of  the  rate  constant  to  a  rate  constant  for  xenon,  Olander 
states  that  in  the  case  of  an  atomic  species  other  than  a  vacancy,  all  one  needs  to  do 
is  replace  the  vacancy  diffusion  coefficient  in  the  above  equation  with  the 
appropriate  atomic  diffusion  coefficient  Interestingly,  had  the  appropriate  atomic 
diffusion  coefficient  been  used  in  dander’s  development  of  the  above  equation  the 
vacancy  fraction,  would  have  appeared  in  the  denominator  of  the  rate  constant 
This  happens  because  of  dander’s  definition  of  the  diffusion  coefficient  for  an 
atomic  species,  i.e. 

A  -  (  49  ) 

versus  the  definition  for  the  vacancy  diffusion  coefficient,  i.e. 

D^^alw  (50) 

where  w  =  the  jump  frequency  of  the  species.  One  possible  reason  for  this 
discrepancy  is  that  an  atomic  species  interacting  with  another  species  (gas  atom, 
dislocation,  bubble,  grain  boundary)  will  essentially  ‘see’  a  vacancy  fraction  of  one 
due  to  the  vacancies  associated  with  the  secondary  interacting  species.  Certainly, 
the  vacancy  fraction  would  be  expected  to  be  greater  than  in  the  bulk  matrix.  If  a 
vacancy  fraction  of  one  is  not  the  case,  then  assuming  it  to  be  so  will  give  the 
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maximum  rate  of  interaction.  Making  appropriate  substitutions  and  notation 
changes,  the  rate  equation  becomes 


where  =  diffusion  coefficient  of  single  xenon  atoms  and  z, , 
constant  for  nucleation  is  given  by 


(51) 

=  z^  =  84.  The  rate 


(52) 

leaving  a  nucleation  rate  of 

Nucleation  Rate  -  (  53  ) 

where  was  multiplied  by  two  to  account  for  both  gas  atoms  being  mobile  [5]. 

Most  recent  researchers  have  expressed  the  nucleation  rate  from  a  less 
theoretical  approach,  based  on  the  work  of  Chandrasekhar  [34]  who  developed  rate 
equations  for  use  in  colloidal  chemistry.  Rest  [17]  expresses  the  nucleation  term  as 
Nucleation  Rate  =  \6nF^D^Cl  (  54  ) 

where  Fj  =  probability  that  two  gas  atoms  which  come  together  actually  stick  and 
R,=radius  of  volume  occupied  by  the  intragranular  gas  atom.  Hayns  and  Wood 
[97]  use  a  similar  expression  that  does  not  include  the  sticking  factor  and  is  lower  by 
exactly  a  factor  of  two. 

An  effective  radius  of  a  gas  atom,  R,,  can  be  calculated  by  assuming  an 
atomic  volume  Q=  djA  and  equating  this  with  the  volume  of  a  sphere,  |  Jr  R?, 
giving 


A 


a  i - =  0. 

“V16jr 


69336a„ 


(55) 


This  gives  a  reaction  rate  (with  F,=  l)  based  on  Chandresakhar’s  work  of 


Nucleation  Rate  »  34.S52\a^D^Cf 


(56) 
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which  is  approximately  20%  less  than  that  of  Olander  (45%  if  a  combinatorial 
number  of  126  is  used).  Esteves  [68]  has  compared  Chandresakhar's  equations  to 
more  exact  (and  complex)  theoretical  rate  equations  [98]  and  concluded  that  both 
are  of  comparable  accuracy.  Given  the  relative  closeness  of  the  two  results,  the 
degree  of  uncertainty  currently  present  in  fuel  modeling,  and  the  conclusions  of  the 
above  study,  it  is  reasonable  to  utilize  Chandresakhar’s  equations  for  bubble 
nucleation.  Since  they  are  also  more  efficient  and  straightforward  in  fiiel  modeling 
codes,  their  use  in  REDSTONE  seems  more  than  justified.  It  will  be  seen  later  that 
this  has  advantages  as  bubble  coalescence  and  atom  capture  can  be  treated  similarly. 
A  sticking  factor  is  not  used  in  program  REDSTONE  as  its  use  by  Rest  appears  to 
be  as  a  fitting  parameter  to  adjust  his  code  to  actual  data 

The  previous  discussion  of  nucleation  was  based  on  theory  involving  random 
effects.  In  a  temperature  gradient,  atoms  and  bubbles  will  experience  an  additional 
force  pushing  them  in  a  single  direction  at  a  speed  proportional  to  the  temperature 
gradient  and  their  size.  For  nucleation,  no  additional  biased  terms  are  reeded  since 
all  gas  atoms  at  a  particular  location  are  pushed  at  the  same  velocity.  Only  random 
effects  will  contribute  to  the  nucleation  of  bubbles 

4.3.3  Coalescence  and  Capture 

Bubble  coalescence  and  gas  atom  capture  are  expressed,  as  for  nucleation, 
through  a  rate  equation.  In  this  case  two  concentrations  must  be  considered,  one  for 
the  bubble  and  one  for  either  the  gas  atom  or  other  bubble,  i.e. 

Coalescence  Rate  =  ^iu,QC„  (  57  ) 

where  C|j,C„  =  concentrations  of  bubbles  with  ‘k’  and  ‘n’  atoms  respectively  and 
k^=  rate  constant  for  the  reaction.  Again,  analysis  of  this  equation  will  be 
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approached  first  from  theoretical  considerations  and  then  from  the  equations  of 
Chandrasekhar,  the  form  typically  used. 

Theoretical  analysis  of  the  rate  constant  has  been  performed  by  Olander  [5]. 
The  method  of  analysis  he  uses  depends  on  the  predominant  process  controlling  the 
reaction.  Coalescence  of  reacting  species  (bubbles,  atoms,  or  vacancies)  which  are 
of  comparable  size,  and  thus  comparable  mobilities,  is  considered  to  be  reaction  rate 
controlled,  as  was  the  nucleation  rate  developed  in  the  previous  section.  When  one 
of  the  species  is  much  larger  than  the  other,  the  process  is  diffusion  limited  due  to 
the  establishment  of  a  concentration  gradient  of  the  more  mobile  species  around  the 
larger  species.  Between  these  two  extremes  a  combination  of  these  two  processes 
exists.  Interactions  between  gas  atoms  and  bubbles  will  be  treated  as  reaction  rate 
controlled  as  bubble  sizes  are  expected  to  be  small  in  the  scenario  here.  Interactions 
between  bubbles  will  be  treated  likewise.  Interactions  between  bubbles  and 
vacancies/interstitials  are  treated  as  diffusion  controlled,  both  theoretically  and  in 
practice,  as  vacancies  and  interstitials  will  have  much  higher  diffusivities  than  the 
bubbles  involved  in  the  interaction.  Diffusion  controlled  reactions  will  be  discussed 
in  the  section  dealing  with  vacancy- bubble  interaction. 

In  general,  bubble  interaction  rates  based  on  reaction  rate  control  are 
cumbersome  to  calculate  from  purely  theoretical  concerns.  Partly,  this  is  because 
one  must  calculate  analogous  terms  to  the  jump  distance  and  jump  frequency  used 
for  atomic  species.  The  difficulty  in  calculating  these  terms  involves  the  more 
general  question  of  what  constitutes  the  formation  of  the  coalesced  bubble. 

Although  the  time  required  for  a  single  coalescence  is  usually  considered  to  be 
negligible,  it  is  not  instantaneous  and  involves  the  movement  of  many  atoms  from 
start  to  finish.  This  complication  manifests  itself  in  the  calculation  of  the 
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combinatorial  number,  used  in  the  rate  constant.  However,  the  interaction  of  a 
single  gas  atom  (treated  as  an  atomic  species)  and  a  stationary  bubble  (gas  atom 
capture)  is  simple  enough  that  a  combinatorial  number  can  be  calculated  using  a 
modification  of  Olander’s  analysis  of  vacancy- vacancy  interactions.  In  the 
following  analysis  gas  atom  capture  by  a  bubble  will  be  considered,  with  the 
modification  that  the  vacancy  fraction  will  be  taken  as  one. 

The  rate  constant  for  gas  atom  capture,  analogous  to  that  for  nucleation,  is 
given  by 


(  58) 


where  Q  =  atomic  volume,  =  lattice  constant,  z,„  =  combinatorial  number,  and 
£)j=  diffusion  coefficient  of  the  single  atom.  Because  only  D  ,  is  present  in  the 
above  equation,  it  is  implicitly  assumed  that  the  bubble  is  stationary.  In  order  to 
calculate  a  rate  constant  for  two  mobile  species,  it  is  necessary  to  calculate  A:„j  and 
add  this  to  It  is  for  this  reason  that  the  nucleation  rate  constant  in  the  previous 
section  was  multiplied  by  a  factor  of  two  (yt,j+A,j).  As  discussed  previously,  the 
calculation  of  is  complex  as  it  requires  the  probability  of  a  bubble  ‘jumping’  next 
to  a  gas  atom.  The  bubble  has  been  assumed  stationary  for  this  reason.  However, 
by  considering  only  an  acceptable  value  for  comparison  can  still  be  obtained 
since  the  magnitude  of  will  be  smaller  than  due  to  the  gas  atom  being  the 
more  mobile  species. 

The  development  of  the  rate  constant,  following  Olander,  begins  by  first 
considering  the  rate  of  formation  of  bubbles  containing  ‘n+1  ’  atoms  through  capture 
of  a  single  gas  atom,  i.e. 

Interaction  Rate  =  /Ji+i^i 


(  59) 
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where  =  probability  that  a  gas  atom  jumps  into  a  site  nearest  neighbor  to  the 
bubble.  A  gas  atom  that  jumps  into  a  nearest  neighbor  site  is  assumed  to  become 
part  of  the  bubble.  The  probability  P^^  depends  on  the  number  of  nearest  neighbors 
as  well  as  the  probability  per  second  that  a  gas  atom  jumps  into  a  particular  site. 

The  number  of  nearest  neighbors  to  a  bubble  has  previously  been  calculated  to  be 
Nearest  Neighbors  =  8. 40 128 .  (  60  ) 

Thus,  the  probability  per  second  that  a  gas  atom  jumps  into  a  nearest  neighbor  site  is 
given  by 

/i„-8.40128;lR.V/>,  (61) 

where  P^  =  probability  per  second  that  a  gas  atom  jumps  into  a  particular  site.  The 
probability  P^  can  be  obtained  by  considering  the  number  of  sites  next  to  a  specific 
nearest  neighbor  site  which  are  not  themselves  nearest  neighbors.  One  must  also 
consider  the  probability  that  a  particular  site  is  occupied  by  a  gas  atom  and  the  jump 
frequency  of  the  gas  atom  in  a  particular  direction. 

The  number  of  sites  next  to  a  nearest  neighbor  will  range  between  four  and 
seven.  Four  sites  are  present  when  considering  a  planar  bubble  surface  (as  for  a 
large  bubble),  while  seven  sites  are  present  in  the  case  of  nucleation.  The 
probability  that  a  particular  site  is  occupied  by  a  gas  atom  is  given  by  C,i3.  The 
jump  frequency  of  gas  atoms  in  the  matrix  can  be  obtained  from  theory  for  the 
diffusion  coefficient  =  1),  i.e. 

A  =  a>  (62) 

where  w  =  jump  frequency  of  xenon  atoms.  Rearranging,  one  gets  the  jump 
frequency  in  terms  of  measurable  quantities. 


(63) 
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Putting  these  quantities  together  gives 
P„^i-[8.40128;r^](4C,fi)[^]  =33.6051jr-^Z),C,r2  . 


(64) 


The  rate  of  gas  atom  capture  for  bubbles  of  ‘n’  atoms  is  then 
Capture  Rate  -  33. 6Q5\k^ I\C^QC^ 
giving  a  rate  constant  of 

-  £ls^  „  33. 605  Iji^DyQ 
and  a  combinatorial  number  of 


(65) 


(66) 


z,„-33.6051jr^ 


(67) 


The  value  for  for  UN  is  approximately  4.889x10’*  cm.  Substituting  this  into  the 
rate  equation  gives 

Capture  Rate  -  1.7184  x  \Q^ nP^D^C^C„  .  (  68  ) 

As  with  nucleation,  many  investigators  [17,97]  handle  random  gas  atom 
capture  following  Chandrasekhar  [34],  i.e. 

Rate  -  A7t(l\  +  £)„)(R,  +  R„)C,C„  .  (  69  ) 

The  form  of  this  equation  is  identical  to  that  used  for  bubble  coalescence,  and 
similar  to  that  for  bubble  nucleation.  It  thus  provides  a  uniform  theoretical  treatment 
of  gas  atom  and  bubble  interactions.  Both  bubble  and  gas  atom  motion  are  included, 
i.^  neither  species  is  assumed  stationary.  If  the  bubble  is  assumed  stationary 
(Z)„=0)  the  rate  equation  becomes 

Rate  -  4;r£>,(R,  +  /?„)C,C„  .  (  70  ) 
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This  equation  differs  from  that  derived  using  Olander’s  approach  through  a  linear 
versus  quadratic  dependence  on  bubble  radius.  An  inspection  of  the  equation  based 
on  Chandrasekhar  reveals  that  as  the  bubble  becomes  much  larger  than  a  single 
atom,  the  theoretical  form  of  the  diffusion  limited  reaction  is  achieved,  i.e. 

Diffusion  Limited  Rate  =■  .  (  71  ) 

The  equations  of  Chandrasekhar,  at  first  glance,  appear  vastly  different  from 
the  reaction  limited  equations  based  on  a  purely  theoretical  treatment  (with  its 
associated  assumptions).  A  more  direct  comparison  can  be  made  by  first 
considering  nucleation  and  substituting  for  R,  and  /?„,  which  are  both  approximately 
3.072x  lO"*  cm.  In  this  case,  the  theoretical  equation  (with  Z)„  assumed  to  be  zero) 
becomes 

Nucleation  Rate  -  1.62  x  10"’;r£\C,C„  (  72  ) 

while  the  equation  according  to  Chandrasekhar  becomes 

Nucleation  Rate  -  2.46  x  10"’;eD,C,C„  .  (  73  ) 

The  two  rates  are  comparable  and  certainly  within  any  error  associated  with  the 
assumptions  made.  If,  instead,  a  bubble  of  25x10-*  cm  radius  is  assumed  for 
then  the  rates  become 

Theoretical  Rate  =  1.07  x  10"^;rD,C,C„  (  74  ) 

and 

Chandresakhar  Rate  =  1. 12  x  10~^;r£),C,Cn  .  (  75  ) 

The  difference  between  rates  is  now  an  order  of  magnitude  and  increasing  as  bubble 
size  increases.  Another  shortcoming  of  the  theoretical  analysis  was  that  the  bubble 
was  assumed  to  be  stationary,  at  least  relative  to  the  gas  atom,  which  is  certainly  not 
the  case  for  analysis  just  presented.  Thus  the  accuracy  of  the  theoretical  treatment 
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as  applied  is  questionable,  even  for  the  regime  where  reaction  rate  control  would  be 
dominate.  It  should  be  apparent  from  the  above  analysis  that  the  theoretical 
approaches  for  bubble  and  gas  atom  interactions  are  not  applicable  over  the  full 
range  of  bubble  sizes.  The  approach  using  the  equations  of  Chandrasekhar  does 
seem  to  accommodate  the  full  range,  with  acceptable  agreement  for  reaction  rate 
control  and  convergence  to  the  theoretical  limit  for  diffusion  control.  The 
conclusion  can  be  made  that  the  theoretical  approach  as  applied  by  Olander  is 
inaccurate  for  the  full  range  of  bubble  sizes  and  too  difficult  to  use.  The  approach 
taken  from  the  work  of  Chandrasekhar  provides  a  viable  alternative  to  describe 
interactions  over  the  entire  range  of  bubble  sizes.  The  reaction  rate  for  interaction  of 
two  bubbles,  according  to  Chandrasekhar,  would  be 

Rate  -  (  76  ) 

which  is  also  the  same  form  resorted  to  by  Olander  in  applying  coalescence  to  his 
own  program  and  is  the  form  chosen  for  use  in  REDSTONE.  For  self-coalescence 
loss,  this  rate  is  multiplied  by  two  since  two  species  of  ‘n’  size  are  removed  for  each 
reaction. 

For  biased  effects  due  to  a  temperature  gradient,  additional  terms  are 
required  for  any  two  species  of  different  size.  Bubbles  of  the  same  size  would  not 
interact  for  the  same  reason  there  is  no  biased  nucleation.  The  force  experienced  by 
a  bubble  in  a  temperature  gradient  is  due  to  the  difference  in  surface  atom  energies 
between  the  ‘hot’  and  ‘cold’  sides  of  the  bubble.  The  atoms  on  the  hot  side  will 
have  more  energy  allowing  more  atoms  to  escape  from  the  surface  zmd  travel  along 
the  inside  of  the  bubble.  Some  of  these  atoms  will  come  to  rest  on  the  cold  side 
where  they  will  lose  their  energy  and  remain.  The  result  will  be  a  net  flux  of  surface 
atoms  to  the  cold  side  of  the  bubble  and  a  migration  of  the  bubble  up  the 
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temperature  gradient.  Biased  bubble  motion  is  treated  as  if  a  force  is  acting  on  the 
bubble  as  a  whole,  and  is  proportional  to  the  bubble  size  and  the  temperature 
gradient  within  the  bubble.  Because  of  the  lower  thermal  conductivity  of  the  gas 
within  the  bubble,  the  temperature  gradient  within  the  bubble  is  higher  than  in  the 
bulk  matrix.  For  a  spherical  bubble,  the  effective  temperature  gradient  on  the 
bubble,  is  related  to  the  temperature  gradient  in  the  bulk  matrix,  Vr,  by 

vrb-|vr.  (77) 

The  velocity  of  a  bubble  in  a  temperature  gradient  is  given  by  the  Nemst-Einstein 


(78) 


where  k  =  Boltzman  constant,  Z),  =  Surface  diffusion  coefficient,  and  the  force,  F^,  is 
given  by 


(79) 

\  <  )  T 

where  is  the  activation  energy  for  surface  atoms  breaking  loose  from  the  surface. 
Although  £)g  has  been  measured  for  UN,  ^  has  not.  However,  Olander  [5]  notes 
that  the  value  recommended  by  Maiya  [79]  for  UO^  is  approximately  80%  of  the 
heat  of  vaporization.  This  result  is  reasonable  considering  that  jumps  of  surface 
atoms  involve  atoms  which  almost  leave  the  surface  permanently.  Thus,  an  estimate 
of  QI  for  UN  can  be  obtained  by  assuming  it  to  be  80%  the  heat  of  vaporization  . 

For  UN,  this  value  has  been  reported  as  1. 1250x  lO"'*  ergs  [99],  corresponding  to  a 
^  of  9.0003X  lO"*^  ergs. 
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In  a  given  temperature  gradient,  it  is  the  size  of  the  bubble  that  determines 
the  velocity.  The  coalescence  rate  will  therefore  depend  directly  on  the  relative 
velocities  of  the  two  bubbles,  as  well  as  the  cross  section  for  interaction,  i.e. 

Biased  Coalescence  Rate  -  .  (  80  ) 

4.3.4  Resolution 

The  rate  of  resolution  of  gas  atoms  from  bubbles  is  commonly  treated  by 
applying  a  resolution  parameter,  b,  to  the  bubble  concentration,  i.e. 

Resolution  Rate  =•  bCj^N^^  (81) 

where  C„  =  concentration  and  =  number  of  gas  atoms  in  the  bubble. 

Experimental  measurements  [100,101]  have  placed  the  resolution  parameter  for  UO-, 
in  the  range  1.8x10^  <  |h|  <  3.6x10"*.  Theoretically,  the  resolution  parameter  has 
been  described  by  a  variety  of  different  methods  based  on  various  theories.  One 
theory  of  resolution  [102]  attributes  the  process  to  the  production  of  a  thermal  spike 
from  a  passing  fission  fragment  in  the  vicinity  of  a  gas  bubble.  The  thermal  spike 
creates  temporary  disorder  in  the  lattice  surrounding  the  bubble  completely  mixing 
fission  gas  atoms  with  matrix  atoms.  The  simplest  model  of  this  theory  [103] 
proposes  that  all  bubbles  intersected  by  a  fission  fragment  track  are  entirely 
destroyed.  Obviously,  this  model  works  best  for  small  bubbles.  The  resolution 
parameter  in  this  case  is  given  by 

b  -  InP^ix^F  (  82  ) 

where  /?„  =  radius  of  bubble  with  ‘n’  atoms,  =  distance  travelled  by  a  fission 
fragment  during  slowing  down  from  birth  energy,  and  =  the  bubble  cross 
section  for  annihilation.  Because  whole  bubbles  are  destroyed  in  the  above  theory, 
the  loss  rate  of  bubbles  due  to  resolution  is  given  by  hC„  rather  than  hC^V„. 
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A  modification  of  the  above  model  [104]  suggests  that  the  cross  section  for 
aimihilation  is  not  the  bubble  cross  section,  but  a  cross  section  determined  by  a 
cylinder  of  disorder  created  along  the  fission  fragment  track.  The  cross  section  for 
annihilation  would  then  become  -  Rf),  where  Rf  is  the  radius  of  the  damage 
cylinder.  Complete  annihilation  of  a  bubble  would  occur  if  it  was  totally  within  this 
cylinder.  For  bubbles  larger  than  this  cross  section,  a  passing  fission  fragment 
would  simply  decrease  the  size  of  the  bubble. 

Another  proposed  theory  of  resolution  [105]  is  that  a  passing  fission 
fragment  blasts  off  chunks  of  the  matrix  onto  the  opposite  wall  of  a  bubble  through  a 
pressure  pulse,  trapping  fission  gas  in  the  process.  Full  conditions  for  the  resolution 
processes  above  have  been  discussed  by  Blank  and  Matzke  [8].  However,  the  above 
theories  probably  do  not  describe  resolution  in  a  fuel  like  UN,  which  has  a  higher 
thermal  and  electrical  conductivity  than  UOj  [8,106].  For  UN  the  resolution 
constant  is  expected  to  be  at  least  an  order  of  magnitude  less. 

Other  resolution  models  [10,107,108],  based  on  theory  attributed  to  Nelson, 
propose  a  knock-on  process  whereby  the  fission  fragment  physically  knocks 
individual  gas  atoms  out  of  the  gas  bubble  into  the  matrix.  In  this  theory,  the 
resolution  parameter  is  given  by  [5] 

^min 

b-  j4,{E„)a{E„.T)dE^.2o{Er.T^. 

e,r 

where  Ef(  =  the  energy  of  a  fission  fragment  at  some  point  in  the  slowing  down 
process,  =  the  maximum  energy  of  a  fission  fragment,  =  the  energy 
spectrum  of  the  fission  fragment  flux,  =  minimum  energy  required  to 
completely  redissolve  the  gas  atom,  o(Eff,T)  =  transfer  cross  section  for 


(83) 


62 


energy  7lo  T+dT  from  a  fragment  of  energy  This  resolution  parameter  is  highly 

dependent  on  the  minimum  energy  needed  for  resolution.  This  minimum  energy 
will  increase  with  increasing  temperature  since  gas  atom  mobility  at  higher 
temperatures  will  make  it  easier  for  an  escaped  gas  atom  to  diffuse  back  to  the 
bubble.  Olander  calculates  the  above  resolution  parameter  using  values  of 
£^**=67MeV,  fi  {{=6  micrometers,  and  two  values  of  300  eV  and  1  keV.  The 
calculated  values  are  A=l.  IxlO'^^Fand  Z>=4.1xlO'‘’Ffor  Tj^j^=300  eV  and  1  keV, 
respectively.  These  calculated  values  are  about  an  order  of  magnitude  lower  than 
experimental  values  in  UO2.  Nelson  [10]  tried  to  correct  for  this  by  estimating  the 
additional  effect  of  resolution  from  secondary  knock-ons.  However,  in  a  fuel  such 
as  UN,  it  is  not  clear  that  secondary  knock-ons  would  be  produced  as  easily  as  in 
UO2.  For  UC,  the  theoretical  single  knock-on  resolution  rate  has  been  shown  to  be  a 
satisfactory  estimate  of  the  experimental  resolution  rate  [9].  Since  UN  has  similar 
electrical  properties,  the  same  conclusion  is  made  for  UN.  Thus,  values  of  the 
resolution  parameter  used  in  REDSTONE  come  from  the  theoretical  calculated 
values  of  Olander.  Resolution  at  low  temperatures  (<500  K)  utilizes  the 
rn^=300eV  value,  while  resolution  at  high  temperatures  (>  1 500  K)  utilizes  the 
r„,i„=l  keV  value,  with  prorated  values  between  these  two  extremes. 

The  above  resolution  parameter  is  all  that  is  required  for  small  bubbles.  For 
large  bubbles,  the  energy  imparted  to  a  gas  atom  in  the  center  of  the  bubble  may  not 
be  enough  to  force  the  atom  past  other  gas  atoms  and  out  of  the  bubble.  Only  gas 
atoms  located  within  a  distance,  d,  of  the  surface  will  have  enough  energy  to  be 
ejected  from  the  bubble.  Nelson  estimates  this  distance  to  be  about  15  A,  increasing 
with  decreasing  bubble  gas  density  according  to  the  relationship. 
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(84) 


where  B  =  the  van  der  Waals  constant  for  xenon  and  =  the  gas  bubble  density. 
From  the  van  der  Waals  equation  of  state  and  the  proportion  of  the  outer  shell  of 
thickness,  d,  to  the  total  bubble  volume,  a  resolution  efficiency  for  large  bubbles  is 
obtained,  i.e. 

-.3 


1-15 


J_ 

K 


(85) 


where  =  the  bubble  radius  ( A),  5  =  the  van  der  Waals  constant,  k  =  Boltzman 
constant,  and  y=  surface  energy.  Thus  the  resolution  rate  for  bubbles  with  /?n>15  A 
is 

Resolution  Rate  =  .  (  86  ) 

Since  resolution  introduces  single  gas  atoms  back  into  the  matrix,  the  summation  of 
all  resolutions  of  bubbles  provides  a  source  term  for  single  gas  atoms. 


4.3.5  Diffusion  Loss  Due  to  a  Temperature  Gradient 

When  the  grain  is  subjected  to  a  temperature  gradient,  gas  atoms  and  bubbles 
will  be  forced  in  a  single  direction  without  regard  to  the  symmetry  of  the  gram.  For 
long  term  operation,  this  introduces  serious  complications  to  the  modeling  of  the 
fuel  grain,  i.e.  one  must  develop  additional  equations  and  methods  to  handle  the  3-D 
problem.  Because  of  the  difficulty  of  modeling  the  3-D  equations,  unique  solving 
methods  are  required.  Fortunately,  the  scenario  to  be  modeled  here  introduces  few 
biased  effects  as  a  result  of  1)  the  essentially  negligible  temperature  gradient  during 
long  term  operation,  and  2)  the  short  duration  of  operation  when  a  temperature 
gradient  becomes  important.  The  simplifications  to  the  3-D  case  to  be  made  here 
are  designed  to  retain  as  much  of  the  biased  effects  as  possible,  but  still  maintain 
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1-D  symmetry  in  the  grain.  Biased  effects  have  already  been  introduced  into  the 
coalescence  terms.  In  this  section,  a  method  will  be  described  for  accounting  for  gas 
and  bubble  loss  to  the  grain  boundary  due  to  the  temperature  gradient. 

In  relation  to  a  spherical  fuel  grain,  the  effect  of  a  biased  force  will  be  to 
shift  all  bubbles  of  the  same  size  a  distance  ‘/i’  of  a  time  period  At.  This  distance  is 
determined  by  the  velocity  of  the  bubbles.  In  the  case  of  different  sized  bubbles,  a 
separate  ‘/i’  would  be  associated  with  each  size.  Other  investigators  [15, 109, 1 10] 
had  tackled  this  asymmetry  problem  by  simply  calculating  the  amount  of  overlap 
between  two  identical  spheres  separated  by  a  distance  ‘A’.  The  fraction  of  this 
volume  to  the  total  volume  of  one  of  the  spheres  was  taken  to  be  the  fraction  of 


fission  gas  retained  in  the  fuel  grain.  This  appears  to  be  the  only  method  used  by 
1-D  codes  to  model  transient  release  of  fission  gas  from  the  grain.  The  volume  of 
the  overlapping  regions  is  given  by  [1 1 1] 

where  =  grain  radius.  Dividing  this  by  the  volume  of  the  sphere,  |rr/^,  and 
subtracting  from  one  gives  the  fraction  released,  i.e. 


Release  Fraction  =  1  - 


16Rj 


(88) 


The  distance  ‘/i’  is  determined  by  multiplying  the  bubble  velocity  by  the  time  step. 


The  bubble  velocity  used  is  the  average  across  the  fuel  grain  for  a  particular  bubble 


size.  A  sphere  the  size  of  the  fuel  grain  is  then  shifted  to  determine  the  percentage 
loss  of  bubbles  of  that  size.  The  percentage  loss  is  applied  to  all  nodes  in  the  grain. 
Each  bubble  size  is  considered  separately,  but  elements  and  nodes  are  not.  This 
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relatively  simplistic  approach  was  used  because  of  the  difficulty  of  accounting  for 
asymmetric  movement  of  bubbles  from  element  to  element  and  node  to  node. 


4.4  Bulk  Vacancy  and  Interstitial  Concentrations 

The  processes  governing  bulk  vacancy  change  in  the  fuel  grain  are  described 
by 


dt 


—  aCyCj  - 


k-l 


Cy-(X  exp 


-—(p  -  —  -P 


(89) 


where  the  terms  are  in  order  of  appearance:  diffusion,  source,  loss  to  dislocations, 
and  loss  to  bubbles.  The  corresponding  equation  for  interstitials  is 
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^  -  oyq  +  jCi  -  aQCi  -  Z|  AftC,  -  2 


N 

I 

k-l 


(90) 


where  each  of  the  terms  are  analogous  to  those  in  the  equation  for  vacancies. 
Effects  from  grain  boundaries  are  accounted  for  through  appropriate  boundary 
conditions,  not  as  explicit  terms  in  the  equations.  Each  of  the  above  terms  will  be 
discussed  in  the  following  sections. 


4.4. 1  Production  by  Fission 

The  production  of  vacancies  and  interstitials  via  fission  has  not  been 
experimentally  investigated  for  UN.  For  UOj,  researchers  have  estimated  a 
production  rate  of  10'^Fi3to  5x10^  [1 12,  113]  with  better  agreement  to 

experiment  being  obtained  for  values  close  to  5x10^  FQ.  Since  production  of 
fission  damage  in  UN  is  expected  to  be  lower  than  in  UOj,  one  would  expect  that 
production  of  vacancies  and  interstitials  would  also  be  lower.  An  estimate  of  the 
production  rate  in  UN  can  be  obtained  by  comparing  experimental  values  for 
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fission-enhanced  diffusion  in  both  UN  and  UOj-  Since  fission-enhanced  diffusion  is 
attributed  to  the  vacancies  produced  from  fission,  such  a  comparison  should  provide 
a  basis  to  estimate  vacancy  production  in  UN  from  the  corresponding  value  for  UO,. 
Fission-enhanced  diffusion  was  discussed  in  section  2.3. 2.2.  In  that  section,  it  was 
seen  that  the  diffusion  rate  in  UO2  was  one  to  two  orders  of  magnitude  larger  than 
that  in  UN.  Assuming  the  vacancy  and  interstitial  production  rates  in  UN  are  the 
same  order  of  magnitude  lower,  the  production  rates  for  UN  can  be  estimated  to  be 
=  5x\0^ FQ  .  (91) 


4.4.2  Recombination 

Olander  [5]  has  given  a  detailed  discussion  of  the  interactions  between  point 
defects,  expressing  rate  constants  in  terms  of  combinatorial  numbers.  His  approach, 
discussed  in  previous  sections,  assumed  that  defects  must  jump  into  nearest  neighbor 
positions  to  interact.  However,  Gibson  [114]  has  shown  that  vacancies  have  a  large 
sphere  of  influence  extending  beyond  simple  nearest  neighbors.  Olander  does  not 
specify  what  the  actual  sphere  of  influence  for  recombination  is,  except  to  suggest 
that  the  combinatorial  number  is  approximately  10  times  larger  (approximately  100 
atom  volumes)  than  the  nearest  neighbor  treatment  (12  atom  volumes).  Straalsund 
and  Guthrie  [67]  use  a  sphere  of  influence  ranging  between  12  and  500  atom 
volumes.  Dollins  and  Jurisch  [71]  use  a  combinatorial  number  of  5040,  equivalent 
to  an  atom  volume  of  about  500  atoms.  Using  this  number,  the  loss  rate  due  to 
recombination  is 


Recombination  Rate  =■  aC^Q 


5040i3(Z>.+Z)J^^ 

2 

a. 


(92) 


67 


where  is  the  atomic  jump  distance.  Although  it  is  clear  that  a  simple  nearest 
neighbor  treatment  is  insufficient  for  recombination,  the  exact  sphere  of  influence  to 
use  is  unclear.  For  program  REDSTONE,  the  loss  rate  utilized  by  Dollins  and 
Jurisch  is  adopted  since  it  does  provide  for  a  larger  sphere  of  influence. 


4.4.3  Dislocation  Interaction 

Although  dislocations  are  assumed  to  have  a  negligible  effect  on  bubbles, 
they  are  included  for  vacancies  and  interstitials  because  studies  have  indicated  that 
the  slight  bias  dislocations  have  for  interstitials  provides  the  driving  force  for  void 
growth  [5,67,1 15].  One  of  the  difficulties  in  including  dislocations  is  determining 
their  types  and  densities.  To  complicate  matters,  both  of  these  can  change  over  time 
with  dose  and  temperature.  Theoretical  treatments  have  attempted  to  account  for 
different  dislocation  types  and  applied  them  to  the  behavior  in  structural  materials 
such  as  steel  [1 15,1 16].  However,  dislocation  properties  are  not  well  documented 
for  nuclear  fuels,  especially  under  irradiation  conditions.  Thus,  for  nuclear  fuels,  the 
applied  dislocation  treatment  has  been  much  more  generalized,  with  dose  and 
temperature  effects  not  included.  In  fact.  Rest  [117]  leaves  out  dislocations 
completely.  Typically,  the  interaction  terms  are  described  by 


Vacancy  -  Dislocation  Interaction  Rate  = 


c,-c;exp(^/>,) 


(93) 


and 

Interstitial  -  Dislocation  Interaction  Rate  =  (  94  ) 

where  thermal  emission  of  interstitials  is  considered  negligible  and  Z^,  =  numbers 

characterizing  capture  volumes  associated  with  unit  length,  Pj  =  overall  dislocation 
density,  Q  =  concentrations,  D^,  -  diffusion  coefficients,  and  =  thermal 


(equilibrium)  concentration.  Because  a  good  correlational  predictor  of  dislocation 
types  and  densities  for  nuclear  fuels  is  not  available,  this  same  set  of  equations  is 
used  in  REDSTONE.  The  obvious  limitations  of  this  approach  are  that  the 
dislocation  density  is  held  constant  and  only  one  type  of  dislocation  is  considered. 

The  numbers  and  Zj  are  typically  on  the  order  of  unity.  However,  because 
of  the  larger  strain  field  surrounding  an  interstitial,  Z;  is  greater  than  Z^  by  one  to 
two  percent  This  difference  is  provides  the  driving  force  for  steady-state  void 
growth.  Consistent  with  other  researchers,  the  values  of  Z^  and  Zj  in  REDSTONE 
are  taken  to  be  1.00  and  1.02  respectively. 

The  dislocation  density  was  adopted  from  other  researchers  as  well.  Values 
of  10*  cm/cm^  to  10^^  cm/cm^  have  been  used  [67,1 18,1 19].  Since  UN  is  expected 
to  have  fewer  dislocations  than  other  corresponding  materials  in  an  irradiation 
environment,  the  lower  value  of  10*  cm/cm^  is  used  in  REDSTONE. 

4.4.4  Bubble  Interaction 

Vacancy  and  interstitial  diffusion  rates  are  much  greater  than  corresponding 
bubble  diffusion  rates.  As  a  result,  the  rate  at  which  vacancies  and  interstitials 
interact  with  a  bubble  can  be  calculated  via  diffusion  equations  in  the  vicinity  of  a 
single  stationary  bubble  (diffusion  limited  reaction),  given  by 
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with  the  boundary  conditions 


Q  =  C*  exp 

Ci  -  Cf  -  0 


and 


(97) 


c  =»  c 


Q-ct 


(98) 


plus  initial  conditions,  where  C^C,  =  vacancy  and  interstitial  concentrations, 
£)y,£)j=vacancy  and  interstitial  diffusion  coefficients,  =  vacancy  and 
interstitial  production  rates,  =  capture  volume  associated  with  the  bubble,  and 
R^=  bubble  radius.  The  effects  from  non-equilibrium  bubbles  have  been  included  as 
boundary  conditions  for  vacancies.  A  similar  boundary  equation  would  normally  be 
present  for  interstitials  except  that  the  interstitial  equilibrium  concentration  has  been 
taken  as  zero.  Capture  effects  from  dislocations  are  not  included  in  these  equations 
as  they  are  already  included  in  the  overall  bulk  concentration  equations  and  bubbles 
are  the  dominant  sink  on  vacancies  [93].  Also,  since  it  is  desirable  to  come  up  with 
an  analytic  expression,  recombination  is  not  included  here  so  that  the  equations 
remain  uncoupled.  Recombination  is  also  included  in  the  calculation  of  the  bulk 
vacancy  and  interstitial  concentrations.  Solutions  to  the  resulting  equations  are  still 
not  easily  obtained  in  finite  form  because  of  the  time  derivative.  Additionally, 
establishment  of  the  initial  condition  for  the  time  variable  involves  knowing  the 
exact  profile  around  each  bubble  size  at  the  beginning  of  each  time  step.  This 
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requirement  is  not  practical  since  it  involves  modeling  vacancy  profile  changes  for 
all  bubbles.  A  greatly  simplified  solution  can  be  obtained  by  taking  a 
quasi-steady-state  approach,  i.e.  setting  the  left  hand  side  of  the  above  equations  to 
zero.  Such  an  approach  assumes  that  the  vacancies  and  interstitials  equilibriate 
around  a  bubble  much  faster  than  the  bubble  changes  in  size.  This  assumption  is 
reasonable  since  one  can  argue  that  after  the  initial  establishment  of  bubble  profiles, 
all  bubbles  would  have  vacancy  and  interstitial  profiles  established  around  them. 
Subsequent  coalescences,  resolutions,  etc  would  not  change  these  profiles  much 
from  their  equilibrium  positions.  Re-equilibriation  of  vacancies/interstitials  would 
occur  rapidly  since  only  minor  changes  would  be  required.  In  support  of  this, 
Matthews  [15]  has  concluded  that  a  steady-state  treatment  of  vacancies  is 
appropriate  in  most  situations.  With  this  simplification,  the  equations  become 


0 


dCy\ 

dr\  dr  / 


(99) 


and 


r^  dr  \  dr  / 

with  boundary  conditions  as  before  and  no  initial  conditions. 


(  100) 


Olander  notes  that  the  solutions  to  these  equations  exhibit  two  distinct 
regions  within  the  capture  volume,  each  dominated  by  a  different  process.  In  the 
inner  region  a  diffusion  loss  term  is  dominant  and  production  effects  can  be 
neglected.  In  the  outer  region  diffusion  loss  is  negligible  and  only  production  need 
be  considered.  Thus  approximate  solutions  in  the  inner  region  can  be  found  by 


solving  for  the  equations 
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r^drV  dr) 


(101) 


and 

dr  \  dr  ) 


(102) 


The  outer  region  is  actually  described  by  the  overall  production/loss  equations 
describing  bulk  vacancy/interstitial  concentrations.  The  concentrations  change  little 
in  the  outer  region  and  quickly  attain  the  values  in  the  bulk  matrix. 

Usually,  the  capture  radius,  R^,  is  taken  to  be  much  greater  than  the  bubble 
radius,  R^.  Assuming  this  to  be  true  here,  the  flux  boimdary  conditions  a.tr  =  R^ 
given  earlier  can  be  replaced  with  concentration  boundary  conditions  at  r  =  oo  ,  i.e. 


Q  -  C° 


r  =  Qc 


(103) 


q  -Q 

where  a  ‘b’  superscript  has  been  added  to  denote  concentrations  far  away  from 
bubbles.  These  boundary  conditions  essentially  convert  the  capture  volume  to  an 
infinite  medium.  The  solutions  to  these  greatly  simplified  equations  are. 


(104) 


and 


The  resulting  reaction  rate  for  a  single  bubble  is  given  by 
Reaction  Rate  =  j  J 

where  J  is  the  flux  across  the  bubble  surface  given  by 


(105) 


(  106) 


(107) 
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Making  the  appropriate  substitutions,  the  net  single  bubble  interaction  rates  become 


Vacancy  Interaction  Rate  =  AitR^DJ 


exp 


(108) 


and 

Interstitial  Interaction  Rate  =  (  109  ) 

where  the  ‘b’  superscript  has  been  dropped  from  the  bulk  concentrations.  To  get  the 
interaction  rates  within  a  unit  volume,  one  needs  to  multiply  by  the  concentration  of 
bubbles  in  the  volume,  i.e. 


Vacancy  Interaction  Rate  = 


Q  -  exp 


(110) 


and 

Interstitial  Interaction  Rate  -  AjtR^L\C^C^  .  (Ill) 

Previous  investigators  have  used  the  forms  above  for  use  in  transient  modeling  [71], 


4.5  Bubble  Radius  Adjustment 

A  bubble  vacancy  adjustment  equation  was  described  in  section  4. 1  in  terms 
of  the  various  processes  involved.  In  this  section,  this  equation  will  be  related  to  the 
concentration  equation,  and  then  the  conversions  needed  to  directly  describe  bubble 
radius  changes  will  be  discussed.  Many  of  the  relationships  discussed  here  have 
already  been  discussed  previously. 

4.5. 1  Relation  to  Bubble  Concentration  Equation 

Ignoring  for  the  moment  the  affect  of  the  bulk  vacancy  concentration  on  the 
growth  of  fission  gas  bubbles,  the  resulting  equation  for  bubble  vacancy  growth  is 
related  to  the  bubble  concentration  equations  by  noting 
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RATE  OF  CHANGE  OF 
NUMBER  OF 
VACANCIES 
ASSOCIATED  WITH 
BUBBLE  OF  ‘N’ 
ATOMS 


RATE  OF  GAIN  OF 
BUBBLES  OF  ‘N’  ATOMS 
THROUGH  PROCESSES 
INVOLVING  BUBBLES  OF 
‘K’  ATOMS 

X 

NUMBER  OF  VACANCIES 
PER  BUBBLE  OF  ‘K’ 
ATOMS 

RATE  OF  GAIN  OF 
BUBBLES  OF  ‘N’  ATOMS 
THROUGH  PROCESSES 
INVOLVING  BUBBLES  OF 
‘K’  ATOMS  AND 
BUBBLES  OF  ‘I’  ATOMS 

X 

NUMBER  OF  VACANCIES 
PER  BUBBLE  OF  ‘K’ 
ATOMS 

+ 

NUMBER  OF  VACANCIES 
PER  BUBBLE  OF  ‘  J’ 
ATOMS 

RATE  OF  LOSS  OF 
BUBBLES  OF  ‘N’  ATOMS 

X 

NUMBER  OF  VACANCIES 
PER  BUBBLE  OF  ‘N’ 
ATOMS 

(112) 


Stated  in  words,  this  conversion  is  simply  the  loss  terms  in  the  original  concentration 
equation  multiplied  by  the  number  of  vacancies  per  bubble  for  the  bubble  being  lost, 
plus  the  gain  terms  multiplied  by  the  sum  of  the  appropriate  number  of  vacancies  per 
bubble  for  the  bubbles  being  gained.  To  complete  the  equation,  one  must  add  the 
effects  from  bubble- vacancy  and  bubble- interstitial  interactions  developed  in  section 
4.4.4.  This  equation  describes  the  rate  of  change  of  the  total  vacancy  concentration 
associated  with  bubbles  of  ‘n’  atoms  in  a  unit  volume.  The  total  vacancies 
associated  with  bubbles  of ‘n’  atoms,  is  equal  to  where  C„  is  the  bubble 

concentration  and  is  the  number  of  vacancies  per  bubble  of  ‘n’  atoms.  In 
REDSTONE,  the  bubble  concentration  is  held  constant  during  bubble  volume 
adjustment  calculations.  The  resulting  equation  with  all  terms  included  is 
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where  each  of  the  terms  have  been  discussed  previously. 


4.5.2  Relation  to  Bubble  Vacancy  Equation 

The  bubble  vacancy  rate  equation  presented  in  the  previous  section  is  not  the 
exact  form  used  in  REDSTONE.  The  equation  in  REDSTONE  describes  the  rate  of 
change  of  the  radius  of  a  bubble,  which  involves  a  conversion  of  the  bubble  vacancy 
equation.  The  conversion  is  made  by  noting  that  vacancies  in  a  bubble  are  related  to 
the  radius  by 


(114) 


where  Q  is  the  atomic  volume  for  one  vacancy.  Applying  this  conversion  and 
rearranging  somewhat,  produces  the  following  equation  for  radius  change. 
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In  addition  to  this  conversion,  several  of  the  properties  depend  on  bubble 
radius.  These  include  the  diffusion  coefficients,  production  rates,  bubble  pressures, 
and  resolution  rates.  Each  of  these  dependencies  are  reviewed  in  following  sections. 
These  should  be  considered  in  the  formulation  of  the  finite  element  problem. 
Unfortunately,  inclusion  of  these  adds  even  more  non-linearity  to  the  problem. 
Preliminary  investigations  indicated  that  only  inclusion  of  bubble  pressure 
dependencies  made  a  noticeable  difference  in  the  convergence  behavior. 
Additionally,  calculational  stability  could  be  adversely  affected  by  allowing  the 
terms  used  to  calculate  concentrations  to  change  when  calculating  radii. 

Fortunately,  bubble  pressure  is  only  used  in  the  vacancy  interaction  terms  which 
were  not  used  to  calculate  concentrations.  Disregarding  the  other  radius  dependent 
property  terms  and  keeping  terms  common  with  the  concentration  equation  constant 
during  the  radius  calculation  simplifies  the  finite  element  formation  of  the  radius 
equation  considerably.  The  details  of  the  formulation  will  be  discussed  in  Chapter 


VII. 
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4.5.3  Radius  Dependent  Properties 

Following  is  a  short  review  of  the  radius  dependence  of  fuel  properties  used 
in  program  REDSTONE. 

4.5.3. 1  Diffusion  Coefficient 

The  diffusion  coefficient  has  been  defined  previously  in  section  3. 3. 2. 2  as 


^new^ 


where  the  included  terms  have  been  defined  previously.  The  diffusion  coefficient  is 
seen  to  have  an  dependence  modified  for  larger  bubbles  by  the  sin^()  term  which 
also  depends  on  bubble  radius. 

4. 5. 3. 2  Production  Rate 

The  production  rate  of  bubbles  through  deposition  of  new  fission  gas 
depends  on  the  fraction  of  the  bubble  surface  area  and  the  total  atomic  sites,  i.e. 

A-„  -  (118) 
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where  the  included  terms  have  been  defined  previously.  This  equation  indicates  an 
dependence. 

4. 5. 3.3  Bubble  Pressure 

The  bubble  pressure  is  described  by  the  van  der  Waals  reduced  equation  of 

state, 
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(119) 
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where  the  included  terms  have  been  defined  previously.  The  bubble  pressure  is  seen 
to  have  a  more  complex  radius  dependence.  The  equation  is  still  in  a  form, 
however,  that  can  be  dealt  with  effectively. 

4.5. 3. 4  Resolution  Rate 

For  bubbles  with  radius  greater  than  15  Athe  resolution  parameter  is 
described  by 

/  f  .i\ 

Z,  =6  1- Jl-15  — +  lf— H  (120) 

where  the  included  terms  have  been  defined  previously.  For  smaller  bubbles,  there 
is  no  radius  dependence.  The  resolution  parameter  for  larger  bubbles  is  again  seen 
to  have  a  more  complex  radius  dependence. 
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CHAPTER  V 

CONCENTRATION  FINITE  ELEMENT  FORMULATION 

5.1  Generalities 

Detailed  modeling  of  the  fission  gas  behavior  within  a  fuel  grain  must  take 
into  account  not  only  the  individual  gas  atoms  but  also  all  possible  gas  bubble  sizes 
In  reality,  the  number  of  bubble  sizes  modeled  is  either  capped  at  some  reasonable 
level  or  the  bubbles  are  grouped  into  one  or  more  ‘average’  bubble  sizes  which 
exhibit  properties  which  are  supposedly  an  average  of  the  group.  This  latter 
approach  has  been  the  preferred  method  in  recent  years  as  it  greatly  reduces  storage 
and  running  time  requirements.  Normally,  the  average  group  values  are  calibrated 
to  agree  with  release  results  rather  than  swelling  results,  i.e.  bubble  group  values  are 
based  on  a  bubble  containing  the  average  number  of  gas  atoms  in  the  group  and 
bubble  volume  is  not  conserved.  Programming  the  collapse  of  bubble  sizes  into 
groups  is  only  done  after  multi-bubble  size  calculations  are  performed,  confirming 
the  collapsed  results.  This  has  not  been  done  for  UN.  Since  the  study  here  concerns 
calculating  both  the  release  and  the  swelling,  the  decision  was  made  to  model 
individual  bubble  sizes  rather  than  groups.  Scoping  calculations  had  indicated  that 
the  number  of  different  bubble  sizes  formed  would  be  small  and  therefore 
manageable  within  the  storage  and  time  constraints  of  the  program.  Consideration  of 
the  gas  atoms  as  well  as  ‘N’  bubble  groups  results  in  ‘N+l’  non-linear  equations  that 
must  be  solved  simultaneously.  Because  of  the  complexity  and  non-linearities  of 
these  equations,  numerical  methods  must  be  used  to  solve  them. 
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coefficient  is  constant  here.  This  is  also  the  form  of  the  equation  used  in 
REDSTONE  for  bubble  and  vacancy/interstitial  concentrations.  The  form  used  for 
bubble  volume  adjustment  is  somewhat  different  due  to  higher  non-linearities  and 


will  be  discussed  in  the  next  chapter. 


5.2  Discretization  of  the  Time  Derivative 

Note  that  in  the  general  equation  the  concentration  appears  as  a  function  of 
both  time  and  position.  A  finite  element  formulation  could  be  performed 
incorporating  both  these  variables,  but  this  would  require  nodal  discretization  in 
time  as  well  as  position.  Such  a  discretization  would  require  a  single  solution  of  a 
set  of  equations  describing  all  concentrations  at  all  nodal  positions  and  at  all 
discretized  times.  Such  a  formulation  would  be  unmanageable  for  any  reasonable 
problem.  Instead,  the  time  dependence  is  usually  treated  with  finite  difference 
techniques  and  finite  elements  are  used  for  the  remaining  position  variables.  In  the 
formulation  of  the  current  problem  the  partial  derivative  in  time  is  replaced  with  a 
backward  Euler  approximation  so  that  only  a  finite  element  formulation  with 
position  need  be  done. 

A  backward  Euler  approximation  given  by 
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was  chosen  as  it  does  not  impose  a  maximum  time  step  size  and  thus  is 
unconditionally  stable.  The  notation  used  is  that  the  unknown  quantity  is  at  time 
‘t+At’,  and  the  last  calculated  known  quantity  is  at  time  ‘t’.  The  unconditionally 
stable  feature  is  important  when  dealing  v.'ith  non-linear  problems  because 
divergence  of  time  stepping  schemes  is  not  always  apparent.  In  other  words,  a 
non-linear  problem  using  a  central  difference  or  forward  Euler  time  stepping  scheme 
may  eventually  converge  even  though  the  solution  became  unstable  somewhere  in 
the  process.  As  a  result,  an  answer  can  be  obtained  which  is  totally  meaningless. 
Although  a  backward  Euler  scheme  provides  the  stability  needed  in  a  non-linear 
problem,  it  has  disadvantages  verses  both  the  other  schemes. 

Against  a  central  difference  scheme,  i.e. 


(123) 


which  is  second  order  accurate,  the  backward  Euler  scheme  is  only  first  order 
accurate.  However,  this  trade-off  is  more  than  acceptable  in  ligh  of  the  poor 
accuracy  of  experimental  data  in  UN  which  must  be  input  into  the  calculation.  The 


advantage  the  forward  Euler  scheme. 
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over  the  backward  Euler  scheme  is  not  accuracy,  but  has  to  do  with  the  efficiency  of 
the  resulting  program.  The  forward  Euler  scheme  is  an  explicit  method.  As  a  result, 
for  a  linear  problem  time  stepping  involves  only  updates  to  the  force  vector  The 
use  of  a  forward  Euler  scheme  allows  one  to  calculate  the  stiffness  matrix  only  once 
at  the  beginning  of  the  problem  because  the  non-time  dependent  stiffness  terms 
remain  constant  in  a  linear  problem.  If  a  backward  Euler  is  used,  the  stiffness 
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matrix  must  be  updated  at  every  time  step  due  to  the  dependence  of  this  scheme  on 
However,  in  a  non-linear  problem,  such  is  the  case  here,  there  are  terms  in 
the  stiffness  matrix  which  must  be  updated  at  every  time  step  anyway.  Thus  the 
advantages  of  the  forward  Euler  and  central  difference  schemes  evaporate,  and  given 
the  choice,  the  backward  Euler  is  the  preferred  scheme  to  use  for  non-linear 
problems. 

After  discretizing  the  time  variable,  the  time  derivative  can  simply  be 
incorporated  into  the  already  defined  generic  constants,  i.e.  since  'Q  is  known 
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and 


(126) 
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Thus  in  the  formulation  the  time  derivative  is  incorporated  into  the  generic 
constants,  leaving 

0  -  oyc,  *  a„  +  fl,C,  +  + 1  tiiQ  +  II 
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With  the  time  derivative  discretized,  the  problem  can  now  be  formulated  in  finite 
elements  using  the  Galerkin  method. 


5.3  Development  of  the  Galerkin  Equations 

Transformation  of  the  original  governing  equation  to  the  Galeilcin  equations 
is  presented  here.  In  later  sections,  the  Galerkin  equations  will  be  tailored  to  the 
problem  and  individual  terms  will  be  identified. 
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First  the  general  equation  is  set  equal  to  a  residual  term, 

{R)  -  £>.V"C„  +  a.  +  /3„C,  +  y.C^  +  + 1  .),Q  + 1  |  £,jQCj 
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The  ‘Nq’  Galerkin  residual  equations  (‘Nq’  will  be  defined  later)  are 


dQ  =  0 
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where  is  a  weighting  function  which  will  later  be  identified  as  the  interpolation 
functions.  The  above  provides  for  ‘Nq’  equations  per  bubble  group  (degrees  of 
freedom).  Thus  the  total  number  of  equations  will  turn  out  to  be  ‘Nq’  times  the 
maximum  number  of  bubble  groups.  It  will  be  necessary  to  utilize  integration  by 
parts  to  reduce  the  order  of  differentiation  in  the  above  equation.  This  manipulation 
is  easier  if  the  exact  physical  geometry  is  first  substituted  into  the  integral  equation, 


i.e.  1-D  spherical  coordinates.  Thus, 
dQ  =  4jcr^dr 
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The  4;r  is  a  constant  that  will  appear  in  every  integral  and  thus  will  have  no  effect  on 


the  equations.  It  can  safely  be  dropped.  Applying  integration  by  parts, 
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Note  that  the  second  order  derivative  in  the  original  equation  has  been  reduced  to  a 
first  order  derivative.  The  first  term  on  the  RHS  of  the  equation  is  a  boundary  term 
to  be  evaluated  at  the  boundaries  of  the  element.  For  the  problem  being  looked  at 
here,  the  explicit  boundary  conditions  are  that  the  gradient  in  concentration  is  zero  at 
the  center  of  the  grain  (r  =  0)  and  that  the  concentration  is  specified  on  the  outer 
boundary.  Additionally,  inter-element  boundary  conditions  require  continuity  of 
concentration  and  concentration  gradient.  For  the  inside  element  which  contains  a 
node  at  r  =  0,  the  boundary  term  reduces  to  zero  at  r  =  0  and  thus  requires  only 
evaluation  at  the  outer  boundary  of  the  element.  Upon  assembly  of  the  equations 
and  application  of  the  specified  concentrations,  the  elemental  boundary  terms  cancel 
out  except  where  a  non-zero  concentration  gradient  is  specified  on  a  physical 
boundary  [120].  In  the  problem  developed  here,  there  are  no  specified  concentration 
gradients  and  thus  integration  by  parts  effectively  yields 


with  the  resulting  equations,  which  are  all  equal  to  zero,  becoming 
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5.4  Incrementalization  of  the  Galerkin  Equations 

It  is  desirable  that  the  above  equations  can  be  solved  by  proven  inversion 
matrix  solving  techniques.  To  do  this  the  equation  needs  to  be  put  into  the  form 
‘[K]{C}={F}’  where  [K]  is  the  stiffness  matrix,  {C}  is  the  vector  containing  the 
unknown  values,  and  {F}  is  the  force  or  load  vector.  This  form  implies  that 
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‘[K]  {C }  ’  should  be  a  set  of  linear  equations,  i.e.  [K]  is  constant  with  {C}  being  the 
unknowns  removed  from  the  original  equations.  It  is  clear  that  putting  the  Galerkin 
equations  as  presented  above  into  the  form  ‘[K]{C}={F}’  will  not  result  in  linear 
equations.  However,  several  methods  of  dealing  with  this  situation  are  available. 
Some  of  these  are  reviewed  below.  As  will  be  seen,  these  are  natural  extensions  of 
each  other. 

The  simplest  solution  method  is  to  assign  all  the  constant  and  non-linear 
terms  to  {F}  and  the  coefficients  to  the  linear  teims  to  [K].  The  non-linear  terms 
used  in  {F}  would  be  based  on  values  from  the  previous  iteration,  and  the  stiffness 
matrix  [K]  would  be  constant  during  the  iterations,  i.e. 

[K]  (  136  ) 


where  the  ‘t’  and  ‘t+At’  denote  the  previous  time  and  the  current  time,  respectively 
and  the  ‘i-r  and  ‘i’  denote  the  previous  iteration  and  the  current  iteration, 
respectively.  Because  all  the  non-linearities  are  based  entirely  on  previous  values, 
the  convergence  rate  for  this  method  could  be  too  slow,  depending  on  the  degree  of 
non-linearity.  In  fact,  very  strong  non-linearities  might  prevent  convergence 
altogether.  This  method  does  have  the  advantage  that  [K]  is  constant  and  thus  does 
not  have  to  be  updated  at  every  iteration.  Unfortunately,  in  the  problem  to  be  solved 
here  the  constants  in  [K]  change  anyway  (although  not  at  every  iteration)  due  to 
explicit  calculation  of  some  of  the  coefficients  at  each  time  step. 

Another  method,  which  is  a  variation  of  the  above,  is  to  move  all  the  non- 
linearities  from  {F}  on  the  RHS  to  the  LHS  while  factoring  out  one  C„,  Q,  or  Cj 
(whichever  is  applicable)  to  produce  {C}  and  a  [K]  which  is  partially  dependent  on 
the  unknown  values.  Thus  {F}  would  remain  constant  and  [K]  would  now  have  to 
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be  updated  at  every  iteration  with  values  from  the  previous  iteration.  The  resulting 
set  of  equations  are  denoted  by 

(137) 

As  in  the  first  method,  the  convergence  rate  could  be  too  slow  since  the  driving 
force  towards  convergence  has  simply  been  switched  from  {F}  to  [K], 

Additionally,  [K]  must  now  be  recalculated  at  every  iteration.  Both  of  these 
methods  require  that  the  calculation  of  the  unknowns  at  the  time  of  interest, 
depend  on  the  previously  calculated  values  at  a  previous  time  or  iteration,  qj. 

A  final  method,  and  the  one  chosen  for  the  current  problem,  is  to  utilize  the 
incrementalization  of  C„  and  cause  the  unknown  values  to  be  AC„.  An  incremental 
dependence  is  achieved  by  defining 

AQ  (138) 

where  AC„  =  the  unknown  change  from  the  previous  iteration  and  _  tQ 
Thus  if  traditional  Newton  iteration  is  employed,  the  last  calculated  value  from  the 
previous  time  period  is  all  that  is  needed  to  start  the  iteration.  The  resulting 
equations  are 

For  an  incremental,  non-iterative  method,  the  same  equations  are  used  except  that 
there  is  only  one  iteration  performed  and  thus  the  value  from  the  previous  time 
period,  i.e.  ‘O’th  iteration,  is  always  required.  This  incrementalization  has  the  added 
advantage  that  the  convergence  can  be  easily  observed  since  AC„-*0  as  convergence 
is  achieved.  The  biggest  advantage,  however,  is  that  the  convergence  is  now  driven 
by  changing  terms  in  both  [K]  and  {F}.  This  accelerates  convergence  as  well  as 


86 


provides  additional  stability  for  strong  non-linearities.  In  fact,  this  third 
implementation  can  be  considered  to  be  a  combination  of  the  first  two  methods, 
incorporating  the  best  of  both  (at  least  in  the  problem  to  be  solved  here).  Again  [K] 
must  be  recalculated  at  each  iteration,  but  the  gains  in  stability  and  convergence  rate 
are  considered  to  be  worth  the  increase  in  calculations. 

It  should  be  obvious  from  the  above  discussions  that  the  values  at  time  ‘t+At’ 
are  always  dependent  on  previously  iterated  values  at  time  ‘t+At’  except  for  the  very 
first  iteration,  where  the  values  from  time  ‘t’  should  be  used.  With  this  understood, 
the  writing  of  the  basic  equations  can  be  simplified  by  dropping  the  time  sujjerscript. 
Additionally,  it  is  apparent  that  if  the  unknowns  are  incrementalized  as  defined 
above,  then  only  the  ‘i-l’  designator  appears.  Thus  the  iteration  superscript  can  also 
be  dropped  with  the  understanding  that  no  superscript  refers  to  the  last  known  value, 
i.e.  the  ‘i-l’th  iteration.  Except  for  a  few  obvious  places,  both  these  superscripts 
will  be  dropped  for  further  illustrations,  i.e.  C„  = 

Using  this  modified  notation,  the  incrementalized  concentration  terms  can 
now  be  substituted  into  the  Galeridn  equations.  After  some  rearrangement  the 
Galerkin  equations  become 
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Note  that  the  unknown  quadratic  (non-linear)  terms  have  been  separated  out  on  the 
LHS.  Except  for  these  terms  the  above  equation  would  easily  fit  into  the  form 
‘[K]  {AC}={F}’.  Inspection  of  the  above  equations  will  reveal  that  the  LHS 
controls  the  rate  of  convergence  while  the  RHS  controls  the  accuracy  of  the  result. 
Additionally,  as  convergence  is  approached  the  quadratic  terms  become  negligible 
very  quickly.  Thus  it  can  safely  be  assumed  that  dropping  the  quadratic  terms  will 
only  slow  the  convergence  rate  and  will  have  no  effect  on  the  resulting  solution. 
The  equation  can  further  be  simplified  by  noting  that  the  ‘k’  and  ‘j’  indices  in  the 


double  summation  terms  of  the  LHS  can  be  interchanged  without  affecting  the 
equation.  Making  these  changes  yields  the  equation 


88 


/ 


dr  dr 


AC„ 


ll 

V 

N 

N  N  N 

N  N 

AQ 

k^n 

k«n  k^nj^n 

j#n  k«*D 

lr"dr 


/ 


"  dr  dr  ‘ 

«n  +  finCn  +  PnkQQ  ^  ^kjQC'j 

k«D  k^D  k^nj^n 

(141) 


The  ‘[K]’,  ‘  {AC}’,  and  ‘  {F}’  terms  are  now  readily  identifiable,i.e. 
[K]{AC}  -  {F} 
where 


’"^IK 

AC, 

>  »  . 

> 

-A'nk, 

[aCnJ 

I^nJ 

(142) 
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"dr  dr 


^n+2y„C„+ 


k«D  / 


[r^dr 


(143) 


and 


f 


( 


Vdr 


(144) 


<Pi  PnkQ  +  tlk  +  2  ^kfi  +  2 

J-n  [ 

\  J-n  J-k  / 

The  rest  of  the  matrix  terms  not  identified  are  filled  out  in  a  similar  manner.  The 
terms  above  are  repeated  for  every  weighting  function,  The  final  stiffness  matrix 
[K]  will  consist  of  ‘N^x  N’  rows.  Discretization  of  the  problem  allows  one  to 
determine  the  weighting  functions  and  causes  the  stiffness  matrix  to  become  square 
with  ‘NqX  N’  columns. 
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5.5  Discretization  of  the  Gal^kin  Equations 

So  that  the  Galerkin  equations  can  be  formulated  in  a  form  useful  for  solving 
real  problems,  the  specific  problem  must  be  discretized.  Discretization  is  a 
necessary  prerequisite  to  determining  weighting  functions  and  methods  to  calculate 
values  such  as  ‘Cn’  and  ‘AC„’.  First,  the  problem  is  attacked  in  a  general  form  by 
assuming  physical  elements  with  ‘Nq’  nodes  p>er  element.  Normally,  the  first  and 
last  nodes  (in  a  1-D  problem)  determine  the  physical  boundaries  of  the  element. 
Within  each  element,  the  values  at  the  nodes  will  be  known  or  calculated.  Between 
nodes,  values  will  have  to  be  interpolated  using  weighting  functions 
(akainterpolation  functions),  i.e.  the  ^-.’s  which  were  used  above.  The  requirements 
for  interpolating  functions  are  that  equals  one  at  the  i’th  node  and  equals  zero  at 
all  other  nodes.  There  are  no  requirements  governing  the  behavior  of  the  weighting 
functions  between  nodes,  and  continuity  across  element  boundaries  is  automatically 
taken  care  off  by  the  formation  of  the  Galerkin  equations  themselves.  So  that  the 
nodal  values  can  be  differentiated  from  the  interpolated  values,  a  subscript 
modification  will  be  introduced  which  refers  to  the  local  (elemental)  node  referred 
to,  i.e.  C^j  will  refer  to  a  nodal  concentration  where  ‘n’  indicates  the  bubble  group, 
as  before,  and  ‘i’  indicates  the  node  number.  The  interpolated  values  are  related  to 
the  nodal  values  by 

Q  -  ( 14S  ) 

1  =  1 

and 

Ng 

AC„-|fAC.„  (146) 

1  =  1 


and 
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No 

'•n  =  2^i''n:i  •  (147) 

1  =  1 

These  values  are  substituted  into  the  previously  developed  Galerkin  equations. 
Because  of  the  interpolation  of  AC„,  every  degree  of  freedom  in  an  element 
generates  ‘Nq’  unknowns  to  be  determined  in  the  discretized  Galerkin  equations. 

The  result  is  a  [K]  matrix  that  is  square  of  dimensions  ‘N^x  N’  as  well  as  {AC}  and 
(F)  vectors  of  length  ‘NqX  N’.  Note  that  because  of  the  first  order  gradient  term, 
the  weighting  functions  must  be  continuous  and  once  differentiable.  Had  iteration 
upon  the  original  governing  equation  been  formulated  instead  of  a  Galerkin 
formulation,  the  weighting  functions  would  have  had  to  be  twice  differentiable. 

Once  the  elemental  matrices  and  vectors  have  been  formed,  they  are 
assembled  into  a  global  matrix  and  vectors  which  are  solved  using  common  solving 
routines  such  as  Gaussian  elimination  or  direct  substitution-iteration.  For  the 
problem  here,  Gaussian  elimination  was  chosen  as  the  solving  method. 

5.6  Determination  of  Interpolating  Functions 

Interpolation  functions  are  used  to  interpolate  between  nodes  in  a  problem. 
When  elements  are  defined  within  a  problem,  the  interpolation  functions  need  only 
interpolate  within  an  element.  Values  between  elements  are  accounted  during 
assembly  into  the  global  matrix  through  the  boundary  nodes  connecting  elements. 

In  theory,  one  would  want  interpolation  functions  which  exactly  fit  the  profile  of  the 
unknown  being  modeled.  In  practice,  this  cannot  be  done  since  one  must  know  the 
solution  to  know  the  profile.  However,  a  general  idea  of  the  solution  profile  is 
usually  known,  and  appropriate  interpolation  functions  can  be  chosen.  The  only 
criteria  for  the  i’th  interpolation  function  is  that  it  is  has  a  specified  value  (usually 
one)  at  the  i’th  node  and  is  zero  at  all  other  elemental  nodes.  Other  constraints  can 
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be  imposed,  such  as  requiring  orthogonal  functions,  but  are  not  required.  Since  most 
problem  solutions  are  continuous  and  smoothly  varying,  the  most  often  used 
interpolation  functions  are  based  on  simple  polynomials.  More  complex  functions 
have  been  used  when  the  problem  solution  warrants  it  However,  it  becomes 
increasingly  more  difficult  to  meet  the  nodal  criteria  mentioned  above  with  more 
complex  functions.  Another  method  used  to  keep  the  interpolation  functions  simple 
is  to  keep  the  number  of  elemental  nodes  low.  Thus,  typically,  investigators  use 
only  two  or  three  nodes  per  element  per  dimension.  The  polynomial  interpolation 
functions  for  both  these  sets  are  presented  in  fig.  1 .  The  specific  characteristics  of 
the  problem  solution  here  are  best  modeled  using  three  nodes  per  element  in  the 
outermost  element.  For  this  reason,  program  REDSTONE  was  programmed  to  use 
three  nodes  per  element  for  all  elements. 
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CHAPTER  VI 

BUBBLE  RADII  FORMULATION 


6.1  Generalities 

The  finite  element  formulation  for  the  bubble  radii  equations  involves  more 
detail  than  for  the  concentration  equations  due  to  higher  non-linearities  and  the  fact 
that  some  of  the  fuel  properties  depend  on  bubble  radius.  Specific  dependencies 
were  discussed  in  Chapter  fV  and  included  the  resolution  parameter,  source  term, 
diffusion  coefficient,  pressure,  and  bubble  velocity.  It  was  also  mentioned  in  that 
chapter  that  all  but  one  of  the  dependencies,  pressure,  could  be  left  out  of  the  actual 
finite  element  implementation  in  program  REDSTONE.  The  bulk  of  the 
formulation  method  for  bubble  radii  remains  the  same  as  for  concentrations.  In  fact, 
the  main  difference  is  the  apportionment  of  terms  in  the  governing  equation  to  the 
correct  finite  element  terms.  In  Chapter  IV,  use  was  made  of  a  generic  governing 
equation  includmg  constant,  linear,  and  quadratic  terms.  This  equation  was  then 
incrementalized  before  being  put  into  Galerkin  residual  equations.  Because  of  the 
higher  non-linearities,  the  generic  equation  in  Chapter  IV  is  not  sufficient  for 
describing  the  bubble  radii  equations. 

In  this  chapter,  rather  than  approach  the  formulation  using  a  general 
equation,  each  individual  process  will  be  looked  at  separately.  These  processes  are 
the  same  ones  previously  discussed.  Each  process  will  be  incrementalized,  leaving 
only  constant  and  linear  terms  to  deal  with.  Two  incrementalizations  are  presented, 
one  with  full  radius  dependence  incorporated  and  one  with  the  radius  dependent 
properties  and  concentration  equation  terms  held  constant..  Except  for  pressure,  this 
latter  form  was  the  one  used  in  program  REDSTONE.  Once  incrementalized,  the 
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methods  used  in  Chapter  FV  to  build  the  Galerkin  equations,  the  stiffness  matrix,  and 
force  vector  are  the  same  needed  for  bubble  radii  equations.  The  incrementalized 
bubble  radii  terms  need  only  be  put  into  residual  Galerkin  equations,  with  linear 
terms  being  put  into  the  stiffness  matrix  and  constant  terms  being  put  into  the  force 
vector.  Only  the  diffusion  gradient  term  needs  to  be  specifically  looked  at  in  terms 
of  the  formulation  of  the  residual  Galerkin  equations.  This  occurs  because  the 
bubble  radius,  which  now  modifies  the  diffusion  term,  is  also  a  function  of  position 
in  the  grain. 

Extensive  use  is  made  of  expansion  in  Taylor  series  which  have  been 
truncated  after  first  order  terms,  i.e. 

/(x  + Ax)«/(x)  + .  (148) 

The  purpose  of  using  Taylor  series  is  to  isolate  unknowns  in  linear,  first  order  terms, 
yet  keep  the  non-linear  behavior.  In  practice,  the  formulation  method  used  here 
could  have  been  used  for  the  concentration  formulation  with  exactly  the  same 
results. 
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The  bubble  radii  adjustment  equation  was  presented  in  Chapter  in.  It  is 
reprinted  here  for  reference, 

C„  A'.Q  ^ 

A'„., ^  i 

5  3 

-  ►kljcA-j 

k«n 

k^n 

-  — Z),qc„+— Z),,clQ-qexp[-— (  149) 

^  ^  I  [  ^7’v  JJl 

As  programmed  in  REDSTONE,  the  radii  which  are  changed  during  radii 
calculations  are  the  ones  at  the  very  ends  of  the  terms  and  the  vacancy/interstitial 
interaction  terms.  Other  radii  remain  constant  during  calculations. 

6.2  Incrementalization  of  Individual  Processes 


6.2. 1  Resolution  Loss 

Resolution  loss  is  described  by 
Resolution  Loss  = 

where  the  resolution  parameter,  for  bubbles  with  /q>l  5  A  is  given  by 


br]„  =  b 


3.Yr^ 


where 


(150) 


(151) 


X-1- 


y 

B 


7-15x10*  . 


(152) 


(153) 
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Combining  this  radius  dependency  with  the  loss  term,  incrementalizing  and 
applying  Taylor  series  gives 


Loss.ii„C.A'A+i(.C„A'Jl- 


1^  -iX-Y  T.Xr^ 


(154) 


Note  that  the  use  of  incrementalization  and  Taylor  series  have  converted  the 
resolution  loss  term  into  the  last  known  loss  term  plus  the  correction  term  for  the 
new  time  step.  The  form  above  would  only  be  applied  for  bubbles  with  A 
For  smaller  bubbles,  has  no  radius  dependence.  Such  a  discrepancy  could 
introduce  oscillations  into  the  radius  calculations  if  a  bubble  has  a  radius 
approximately  equal  to  1 5  A  The  corresponding  incrementalized  loss  term  for  small 
bubbles  is 


Loss-ii.,C.A'„(/!„+A;i„)  (155) 

which  is  considerably  simplified.  This  is  also  the  form  used  if  one  leaves  out  the 
radius  dependence  in  the  resolution  parameter  altogether  and  is  the  form  used  in 
program  REDSTONE. 


6.2.2  Resolution  Gain 

Resolution  gain  is  somewhat  more  difficult  to  handle  than  loss  because  of  the 
mixture  of  dependencies  on  two  different  bubble  radii.  However,  the  basic  approach 
for  incrementalization  is  the  same.  The  gain  term  is  given  by 

Resolution  Gain  =  (156) 

K 

where  the  ‘n+l’  subscript  has  been  replaced  by  a  ‘k’  for  simplification.  Using  the 
same  methods  as  for  resolution  loss,  one  arrives  at  (for  /?k^l5  A), 
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Gain  -  +  AC.A^Jfl  -  +  2X^Y^  -  AT^-^Ia/?. 

3  L  L  k^2  k  lc|V  jj^  ^2|  Me 

+  ^QA'JKa^^  -  l)-%  -  2X^Y^  +  2XY- ^  A-l^/^n  ■ 

For  smaller  bubbles  and  retaining  the  radius  dependence  only  partially, 

Gain--jiiC,A'jA-  +  ^A«, 

j  l  »  ^2  “H  ^  “\.J 


which  is  the  form  used  in  REDSTONE. 


(157) 


(158) 


6.2.3  Self-Coalescence  Loss 

Loss  through  coalescence  of  bubbles  with  same  size  bubbles  is  given  by 
Self  -  coalescence  Loss  (  159  ) 

The  diffusion  coefficient  radius  dependence  can  be  expressed  by 

0.-4  (160) 

K 

then  the  loss  term  becomes 

Self  -  coalescence  Loss  (  161  ) 

where 
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force  vector.  Note  also  that  without  the  sine  term,  the  resulting  incrementalized  loss 
term  has  no  dependence  on  bubble  radius. 

If,  instead,  the  diffusion  coefficient  dependence  is  left  out  and  the  original 
term  from  the  concentration  equation  is  held  constant,  then  the  incremental  solution 
becomes 

Self  -  coalescence  Loss  -  -y  +  A/?„]  (  164  ) 

where  only  the  last  term  in  brackets  changes  during  calculation.  This  is  the  form 
used  in  REDSTONE. 


6.2.4  Coalescence  Loss 

Loss  through  coalescence  of  bubbles  with  different  size  bubbles  is  given  by 
Coalescence  Loss  «  +  Rj.)(Z)„  +  Z\)CnQ  .  (  165  ) 

Applying  methods  as  before  with  full  radius  dependence  included  yields. 
Coalescence  Loss  =  +  ^XAi  A)QQ 

l  +  2^jA«i  (165) 

and  with  radius  dependence  partially  left  out. 

Coalescence  Loss  =  +  A/?„]  (  167  ) 

where  only  the  last  term  in  brackets  changes  during  calculation. 
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where  the  ‘n-k’  subscript  has  been  replaced  with  ‘j’  for  convenience.  Applying 
methods  as  before  with  full  radius  dependence  included  yields. 


Coalescence  Gain  =  +  Z\.  +  ^) 

+  ijrCjQ  -4  [•^2'^  -  -  2^0 

+  +  R^+  iRiRj)  A/?k 

+  ^  -^(2^,"  2J^) 

+  + /?^  +  3/?f/2k)]A/?j 

and  with  radius  dependence  partially  left  out. 

Coalescence  Gain  -  j  +  ^)(^j  +  Z\)CjC^ 

where  only  the  last  term  in  brackets  changes  during  calculation. 


(169) 


(170) 


6.2.6  Biased  Coalescence  Loss 

Loss  through  biased  coalescence  with  another  bubble  is  given  by 
Biased  Coalescence  Loss  =  ^jtR^[R^  +  ”  ^klOiQ  (  171  ) 

Bubble  velocity  depends  on  radius  with  the  dependency, 

Vn  -KK-  (  122  ) 

Unfortunately,  incrementalization  of  the  velocity  requires  splitting  the  terms  into 
constant  and  correction  pieces  which  are  then  placed  at  different  locations  in  the 
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stiffness  matrix  and  force  vector.  Keeping  the  absolute  value  incorporated  becomes 
extremely  difficult.  Because  of  this,  the  velocity  terms  are  kept  constant  in  both  the 
radius  dependent  and  radius  independent  formulations.  The  absolute  value  then 
resides  in  the  force  vector  during  calculations.  However,  incorporating  the  other 
radius  dependent  terms  into  the  formulation  gives. 

Biased  Coalescence  Loss  « 

+ Mss;  +  t  /$)!»;  -  rJc.c,A/^ 

+  -  K,|C.QA«,  .  (  173  ) 

Holding  radius  dependent  terms  from  the  concentration  equation  constant  yields  the 
form  used  in  REDSTONE, 

Biased  Coalescence  Loss  «  ~  ) 

where  only  the  last  term  in  brackets  changes  during  calculation. 


6.2.7  Biased  Coalescence  Gain 

Gain  through  biased  coalescence  of  two  smaller  bubbles  is  given  by 


Biased  Coalescence  Gain  «= 


(175) 


where  the  ‘n-k’  subscript  has  again  been  replaced  with  ‘j’  for  convenience.  The 
arguments  for  keeping  the  velocity  terms  constant  in  biased  coalescence  loss  also 
apply  here.  Applying  methods  as  before  to  the  rest  of  the  radius  dependence  in  the 
above  yields 
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Biased  Coalescence  Gain  «« 

+  +  2/4/2f  +  2/?/|^  -  r^lCjQA/Z^ 

+  ^(5^;  +  8^'^  +  3/?f /?,=  +  2R^/i^  V^\Cf^Mi^ 

(176) 

and  holding  the  radius  dependent  terms  carried  over  from  the  concentration  equation 
constant  yields 

Biased  Coalescence  Gain  -  -j  ^  j  j  <J  -  'ipjCt 

X  ^  -,  -2--^  ,-- A/?„+-^A/g,,+— f  A/g.  (177) 

where  only  the  last  term  in  brackets  changes  during  calculation. 


6.2.8  Source  Loss 

Loss  through  deposition  of  another  fission  gas  atom  in  the  bubble  is  given  by 


Source  Loss  *»  j  Ar„Cn/2n  (178) 

where  the  radius  dependence  of  the  source  term  is 

K,-KX-  (179) 

Inclusion  of  this  into  the  incremental  formulation  gives 

Source  Loss  -  j  +  A'„C„A/2„  (  180  ) 

or  if  the  source  radius  dependence  is  left  out. 

Source  Loss  -  j  A„C„[A„  +  A/?„  ]  (181) 


where  only  the  last  term  in  brackets  changes  during  calculation. 
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6.2.9  Source  Gain 

Gain  through  deposition  of  another  fission  gas  atom  in  the  next  smaller 
bubble  is  given  by 


Source  Gain  = 


(182) 


where  the  ‘n-l’  subscript  has  been  replaced  by  ‘k’  for  convenience.  The  radius 
dependence  of  the  source  term  is  as  stated  in  the  last  section.  The  full  radius 
dependent  formulation  is 


Source  Gam 


(183) 


With  the  partially  radius  independent  formulation  being 
Source  Gam  - 

where  only  the  term  in  brackets  changes  during  calculation. 


^  2 ^  AP  I  3 ^  A/? 


(184) 


6.2. 10  Interstitial  Capture 

A  decrease  in  radius  due  to  capture  of  interstitials  is  given  by 


Interstitial  Capture  Decrease  «  . 

K 


(185) 


The  incremental  formulation  of  this  term  is 


Interstitial  Capture  Decrease  - 


1- 


(186) 


This  is  the  form  used  in  REDSTONE.  There  is  not  a  partially  radius  dependent  form 
as  this  term  did  not  appear  in  the  original  concentration  equation. 
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6.2. 1 1  Vacancy  Capture 

An  increase  in  radius  due  to  capture  of  vacancies  is  given  by 


Vacancy  Capture  Increase  =  QD^C^C^ —  . 

K 


(187) 


The  incremental  formulation  of  this  term  is  similar  to  that  for  interstitial  capture,  i.e. 


Vacancy  Capture  Increase  =  (188) 

Again,  there  is  not  a  partially  radius  dependent  form  as  this  term  did  not  appear  in 
the  original  concentration  equation. 


6.2. 12  Vacancy  Emission 

A  decrease  in  radius  due  to  emission  of  vacancies  is  given  by 


1  ^  \  Q  (  2y 

Vacancy  Emission  Decrease  =  QD^C^  — Q  exp)  —  ^  -  P^ 


(189) 


The  incremental  formulation  of  this  term  is  quite  different  from  the  previous  two 
terms  because  of  the  exponential  which  is  dependent  on  both  bubble  radius  and 
bubble  pressure.  The  bubble  pressure  has  a  radius  dependence  given  by 


"  K-y 

where 


(190) 


(191) 


y.lM 

3jt 


(192) 
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After  incrementalization  and  then  Taylor  series  expansion,  one  obtains 

Vacancy  Emission  Decrease  =  Cl  expJ  — [  P  +  —  -  /I  ]  I 

R,  “  K  )\ 


1 

1 

^  1 

1  “ 

-  yfj 

A/2n 


(193) 


which  is  the  form  used  in  REDSTONE. 


Modification  of  the  Diffusion  Term 

The  diffusion  term  for  bubble  radius  adjustment  is  given  by 

Diffusion  Loss  =  r~  . 

3r-  drV  dr  ) 


(194) 


Before  incrementalizing  this  term,  it  will  be  useful  to  first  incorporate  it  into  a 
Galerkin  formulation.  Holding  the  diffusion  coefficient  constant,  the  Galerkin 
residual  equations  for  this  term  are  given  by 

Diffusion  Loss  -  — C .  (  195  ) 

3  J  drV  dr  / 

R 

As  before,  the  4rrcan  be  dropped.  The  above  term  is  different  from  the  formulation 
in  the  previous  chapter  by  the  presence  of  Applying  integration  by  parts  gives 


Diffusion  Loss  « 


^^out 


^  r  2  . 

3  J  dr  dr 


(  196) 


R 

The  first  term  in  brackets  is  evaluated  at  the  element  boundaries.  Like  the 
concentration  gradient,  the  concentration  gradient  times  the  bubble  radius  must  be 
continuous  and  smoothly  varying  across  the  boundary.  There  can  be  no  jumps  in 
bubble  size  at  a  boundary.  As  a  result,  this  term  cancels  out  everywhere  except  on 
physical  boundaries  where  the  gradient  is  specified.  Since  the  gradient  is  never 


» 
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specified  in  the  problem  here,  this  term  can  be  eliminated  entirely.  The  second  term 
involves  the  gradient  of  a  product.  Expanding  this  out  gives  for  the  diffusion  loss 
term 


Diffusion  Loss  =  f f .  (  197  ) 

2  J  dr  dr  2  J  dr  dr 

R  R 

The  result  of  the  R„  has  been  to  add  an  extra  term  to  be  included  into  the 
formulation.  These  terms  can  now  be  incrementalized  and  incorporated  using  the 


techniques  presented  in  Chapter  V. 
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CHAPTER  VII 

PROGRAM  REDSTONE  ORGANIZATION 


7.1  Overview 

Program  REDSTONE  is  organized  as  shown  in  fig.  2.  The  flow  of  program 
control  involves  three  separate  and  one  overall  iteration  cycles.  The  three  separate 
iteration  cycles  involve  calculations  of  bubble  concentrations,  bulk  vacancy  and 
interstitial  concentrations,  and  bubble  radii.  The  overall  iteration  cycle  is  completed 
when  no  additional  iterations  are  required  in  the  bubble  concentration  or  bubble  radii 
steps  during  any  one  overall  iteration.  Time  stepping  is  done  using  time  steps  and 
time  periods.  Periods  define  a  set  of  operational  conditions  for  a  length  of  time. 

Time  steps  are  user  specified  equal  interval  divisions  within  the  time  period  to 
facilitate  calculational  and/or  output  requirements.  Transient  analysis  is 
accomplished  by  stepping  between  the  initial  and  final  period  conditions  based  on 
the  number  of  time  steps  in  a  period.  Within  REDSTONE  is  built-in  the  capability 
to  ‘restart’  from  a  previously  saved  analysis.  This  allows  one  to  set  up  a  specific 
initial  profile  for  subsequent  exploratory  analyses.  A  complete  program  listing  can 
be  found  in  Appendix  1 . 

Input  to  the  program  is  accomplished  through  an  input  file  with  data  in  free 
format  and  identified  by  position  on  a  tine  and  a  unique  line  number.  Comment 
lines  are  allowed  and  identified  by  an  or  ‘C’  at  the  beginning  of  the  line. 
Optionally,  a  restart  input  file  is  also  used  to  input  starting  conditions.  Output  of  the 
program  consists  of  readable  output  at  user-specified  intervals,  and  a  character  based 
file  for  input  to  separate  plotting/data  analysis  schemes.  Additionally,  after  each 
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period,  specific  parameters  are  saved  in  machine  readable  form  for  use  in 
subsequent  restart  analyses  using  program  REDSTONE. 


Fig.  2  REDSTONE  Program  Flow. 
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7.2  Subroutines 

Program  REDSTONE  is  made  up  of  27  subroutines  and  functions.  Each  of 
these  are  listed  below  with  a  short  description  of  their  functions. 

7.2.1  ASSEMB 

Assembles  the  global  matrix  and  force  vector  to  be  solved  for  bubble 
concentrations.  The  global  stiffness  matrix  is  produced  in  a  banded, 
non-symmetrical  form  with  specified  concentration  rows  and  columns  (those  on  the 
grain  surface)  eliminated  but  still  accounted  for.  The  bandwidth  of  the  global  matrix 
IS  twice  the  number  of  bubble  sizes  times  the  number  of  nodes  in  an  element. 


7.2.2  ASSMBR 

Assembles  the  global  matrix  and  force  vector  to  be  solved  for  bubble  radii. 
Again,  the  global  stiffness  matrix  is  produced  in  a  banded,  non-symmetrical  form 
with  specified  radii  rows  and  columns  (those  on  the  grain  surface  and  gas  atoms) 
eliminated  but  still  accounted  for.  The  bandwidth  of  the  global  matrix  is  twice  the 
number  of  bubble  sizes  minus  one  times  the  number  of  nodes  in  an  element.  The 
number  of  bubble  sizes  is  reduced  by  one  because  the  radius  of  single  gas  atoms  is 
not  variable. 

7.2.3  ASSMBV 

Assembles  the  global  matrix  and  force  vector  to  be  solved  for  bulk  vacancy 
and  interstitial  concentrations.  The  global  stiffness  matrix  is  produced  in  a  banded, 
non-symmetrical  form  with  specified  radii  rows  and  columns  (those  on  the  grain 
surface)  eliminated  but  still  accounted  for.  The  bandwidth  of  the  global  matrix  is 
twice  the  number  of  nodes  per  element  times  two. 
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7.2.4  DIFFUS 

Calculates  diffusion  coefficients  for  single  gas  atoms,  bubbles,  vacancies, 
and  interstitials. 


7.2.5  ELKFG 

Assembles  the  elemental  stiffness  matrices  and  force  vectors  of  bubble 
concentrations  for  subsequent  integration  into  the  global  matrix  and  vector  by 
ASSEMB.  The  governing  equations  are  incorporated  in  this  subroutine. 

7.2.6  ELKFR 

Assembles  the  elemental  stiffness  matrices  and  force  vectors  of  bubble  radii 
for  subsequent  integration  into  the  global  matrix  and  vector  by  ASSMBR.  The 
governing  equations  are  incorporated  in  this  subroutine. 

7.2.7  ELKFV 

Assembles  the  elemental  stiffness  matrices  and  force  vectors  of  vacancy  and 
interstitial  concentrations  for  subsequent  integration  into  the  global  matrb  and 
vector  by  ASSMBV.  The  governing  equations  are  incorporated  in  this  subroutine. 

7.2.8  EQURAD 

Calculates  the  equilibrium  radius  of  a  given  bubble  size  with  temperature 
and  gas  atom  content  as  inputs. 

7.2.9  EXPAND 

Adds  the  iterative  solutions  into  the  respective  change  vectors  describing  the 


change  over  a  time  period. 
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7.2. 10  ERROR 

Prints  the  user  specified  error  in  the  output. 

7.2.11  FRTP AG  ^ 

Prints  the  REDSTONE  banner  in  the  output. 

7.2.12  ICOND 

Calculates  arrays  describing  connectivity  between  bubble  groups  at 
elemental  nodes  with  bubble  groups  at  global  nodes.  Used  by  ASSEMB  and 
ASSMBR. 

7.2.13  INPUT 

Assigns  data  read  in  from  the  input  file  to  respective  variables. 

7.2. 14  INTARS 

Calculates  arrays  describing  correlations  of  all  bubble  groups  to  those  which 
are  not  specified.  Used  by  ASSEMB  and  ASSMBR. 

7.2.15  MTXADD 

Performs  matrix  addition  on  arrays. 

7.2.16  OUTPUT 

Prints  results  of  calculations  in  character  format. 

7.2. 17  PROP 

Calculates  fuel,  fission  gas,  and  bubble  properties  at  the  beginning  of  each 
new  time  step.  Also  calculates  initial  bubble  radii. 
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7.2. 18  RADII 

Calculates  new  bubble  radii  based  on  calculation  results. 

7.2.19  RDCARD 

Reads  a  line  of  input  data  from  the  input  file  and  parses  it. 

7.2.20  RESINP 

Reads  the  restart  input  file  when  restart  is  specified. 

7.2.21  RESSAV 

Saves  data  in  machine  readable  form  to  a  restart  file  after  each  time  period. 
Allows  subsequent  runs  to  start  at  a  specific  point  in  time  and  with  an  already 
calculated  concentration  and  radii  profile. 

7.2.22  RNORM 

Calculates  the  norm  of  a  vector. 

7.2.23  SEMBND 

Calculates  the  half-bandwidth  of  the  bubble  concentration  stiffness  matrix. 
Similar  values  for  vacancy/interstitial  concentrations  and  bubble  radii  are  calculated 
in  assembly  routines. 

7.2.24  SETCON 

Adds  the  time  period  changes  into  the  known  values  and  sets  changes  to 

zero. 

7.2.25  SHAPE 

Calculates  finite  element  shape  functions. 
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7.2.26  SOLVE 

Solves  banded,  non-symmetrical  matrices  and  stores  the  result  in  the  force 
vector.  Does  not  check  for  positive  definiteness. 

7.2.27  SWELL 

Calculates  gas  released  and  swelling  using  calculated  results. 

13  Input 

Input  to  program  REDSTONE  is  done  through  80  column  character  records 
residing  in  an  external  file.  An  example  can  be  found  in  Appendix  2.  Each  record  is 
first  scanned  for  or  ‘C’  in  column  one  to  see  if  the  record  is  a  comment.  If  so, 
parsing  is  stopped  and  the  next  line  read  in  for  processing.  If  the  line  is  not  a 
comment  it  is  first  parsed  into  words  (denoted  by  a  breaking  space).  The  first  word 
is  checked  to  see  if  it  is  a  valid  identifying  line  number.  An  invalid  line  number 
causes  program  execution  to  stop  with  an  error.  Numbers  are  used  in  the  first  word 
of  a  line  to  denote  the  type  of  data  present.  The  numbers  used  are  grouped  by 
hundreds  with  subgroups  by  tens.  Thus  a  three  (to  four)  digit  number  is  required 
with  the  first  two  (or  three)  digits  being  significant.  The  last  digit  can  be  used  at  the 
user’s  discretion.  Once  the  data  type  has  been  determined,  a  check  is  performed  to 
insure  the  right  amount  of  data  is  present  on  the  line.  If  so,  the  individual  words  are 
assigned  to  the  respective  variables  and  the  next  line  read  in  for  processing.  After 
processing  all  lines,  the  assigned  data  is  printed  to  the  output  file.  The  specific  types 
of  data  input  are  reviewed  below  by  identifying  line  number. 
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7.3. 1  Card  Number  100 

Contains  user  specified  header  information.  This  information  is  not  parsed, 
but  is  printed  as  is  on  the  output  file.  Only  the  ‘  1’  in  the  line  number  is  significant, 
leaving  numbers  100-199  to  be  used  by  the  user  for  headers. 

7.3.2  Card  Number  200 

Can  contain  a  variety  of  data  depending  on  the  line  number  digit  in  the  tens 
place.  An  identifier  for  whether  British  or  metric  units  are  used  is  identified  by 
‘210’.  The  number  of  Gauss  points  to  use  for  integration  is  identified  by  ‘220’.  A 
0/1  flag  for  whether  the  input  deck  is  reprinted  on  the  output  is  denoted  by  ‘230’.  A 
0/1  flag  to  indicate  whether  a  restart  file  is  to  be  read  is  denoted  by  ‘240’.  Currently, 
the  type  of  units  identifier  is  not  used.  Program  input  and  output  are  in  cgs  units. 
Both  the  ‘2’  and  the  second  digit  in  the  line  number  are  significant. 

7.3.3  Card  Number  300 

Contains  the  number  of  elements,  the  number  of  bubble  sizes  per  node,  the 
number  of  nodes  per  element,  the  number  of  specified  concentrations,  the  number  of 
time  periods,  and  the  maximum  number  of  iterations  to  allow.  Currently,  the 
number  of  nodes  per  element  is  hard  coded  into  much  of  the  program  for  three  nodes 
per  element,  so  ‘3’  should  be  the  only  value  for  this  variable.  The  maximum 
number  of  iterations  is  applied  individually  for  each  iteration  cycle,  not  as  an 
aggregate  for  all  iterations.  Only  the  ‘3’  in  the  line  number  is  significant. 


114 


7.3.4  Card  Number  400 

Previously  used  to  input  the  connectivity  array  of  local  node  numbers  to 
global  node  numbers.  These  are  now  calculated  in  the  program  and  thus  number 
400  is  no  longer  used. 

7.3.5  Card  Number  500 

Contains  the  location  of  nodes  in  centimeters  from  the  center  of  the  fuel 
grain.  A  separate  line  is  used  for  each  node.  Only  the  ‘5’  in  the  line  number  is 
significant. 

7.3.6  Card  Number  600 

Contains  specified  concentration  information.  Each  sjDecified  concentration 
occupies  one  line.  Each  line  contains  a  number  identifying  the  bubble  size,  the  node 
of  the  specified  concentration,  and  the  value  of  the  concentration  in  cgs  units.  Only 
the  ‘6’  in  the  line  number  is  significant. 

7.3.7  Card  Number  700 

Contains  the  radius  of  the  fuel  grain  in  centimeters.  Only  the  ‘7’  in  the  line 
number  is  significant. 

7.3.8  Card  Number  800 

Previously  used  to  input  material  type  identifiers  for  multiple  material  types 
(i.e.,  different  enrichments,  etc).  Currently  much  of  the  program  is  hard  coded  for 
only  one  material  type.  As  a  result  only  one  material  type  is  allowed  and  this  line 
number  is  no  longer  used. 


115 


7.3.9  Card  Number  900 

Contains  the  fuel  porosity,  enrichment,  and  convergence  limit  to  use.  The 
porosity  and  enrichment  are  fractions  less  than  one.  The  porosity  is  not  currently 
used  since  it  is  an  intergranular  property.  Only  the  ‘9’  in  the  line  number  is 
significant. 

7.3.10  Card  Number  1000 

Contains  the  time  length  of  periods  in  seconds,  the  number  of  time  steps  to 
take  in  a  period,  and  the  interval  in  time  steps  which  to  output  results.  One  line  is 
used  for  each  time  period.  Only  the  ‘10’  in  the  line  number  is  significant. 

7.3. 1 1  Card  Number  1100 

Contains  period  operational  data,  one  line  for  each  period.  Each  line 
contains  the  initial  power  density  in  Watts  per  cubic  centimeter,  the  initial 
temperature  in  degrees  Kelvin,  the  initial  temperature  gradient  in  degrees  Kelvin  per 
centimeter,  the  initial  pressure  in  Pascals,  the  final  power  density,  the  final 
temperature,  the  final  tempeiature  gradient,  and  the  final  pressure.  Only  the  ‘  11  ’  in 
the  line  number  is  significant. 

7.4  Output 

Program  REDSTONE  produces  several  forms  of  output.  The  primary  output 
is  a  listing  of  calculational  results  at  user  specified  intervals.  An  example  can  be 
found  in  Appendix  3.  The  first  part  of  this  listing  is  a  ‘REDSTONE’  baimer 
followed  by  the  user  input  header  statements.  This  is  optionally  followed  by  a 
reprint  of  the  input  file.  The  input  data  is  then  summarized  and  printed  in  tabular 
form.  Along  with  the  input  data  are  tables  of  a  few  preliminary  calculations  such  as 
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connectivity  arrays.  All  the  above  information  is  printed  only  once.  The  rest  of  the 
output  listing  consists  of  the  results  of  calculation  at  the  end  of  time  steps  at 
intervals  specified  by  the  user.  Each  of  the  results  listings  begins  with  time  step 
information,  followed  by  the  percent  fission  gas  released  and  percent  swelling.  The 
calculated  bumup  in  several  forms  is  listed  next  followed  by  the  power  density, 
temperature,  and  pressure  used.  Finally,  the  bulk  of  the  results  listing  involves  the 
calculated  bubble  concentrations,  radii,  and  pressures  by  node. 

In  addition  to  the  above  output,  two  other  files  containing  the  calculated 
results  is  written.  One  of  these  consists  of  only  numbers  separated  by  spaces  and  is 
suitable  as  input  to  other  analysis  programs.  The  data  included  in  this  file  are  the 
node  number,  the  radius  of  the  node,  the  bubble  size,  the  bubble  concentration,  the 
bubble  radius,  and  the  bubble  pressure  on  each  line. 

The  second  is  a  machine  readable  restart  file  which  can  only  be  read  by 
program  REDSTONE.  This  file  is  written  automatically  at  the  end  of  each  time 
period  and  overwrites  any  previous  information  in  the  file.  It  contains  the  last 
calculated  time  period  values  only.  This  includes  the  convergence  denominators, 
the  total  fission  gas  atoms  produced,  as  well  as  the  bubble,  vacancy,  and  interstitial 
data,  i.e.  concentrations,  radii,  and  pressures.  As  a  check  for  restart,  this  file  also 
contains  the  number  of  elements,  number  of  bubble  sizes,  and  number  of  nodes  pier 
element  used.  This  information  is  compared  with  the  input  file  upon  restart  to 
ensure  that  a  compatible  problem  is  being  restarted.  A  non-compatibility  results  in 
an  immediate  program  halt  and  an  error. 
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CHAPTER  VIII 
RESULTS 


8.1  Validation  of  Program  REDSTONE  Internals 

Program  results  will  never  be  valid  if  the  program  doesn’t  perform  as 
intended.  So  that  the  internal  consistency  of  REDSTONE  could  be  checked, 
calculated  results  were  compared  with  analytic  solutions  where  possible.  At  least 
two  analytic  solutions  are  available  for  comparison.  Both  are  linear,  single  bubble 
size  solutions.  The  first  is  the  steady-state  solution  and  the  second  is  a  time 
dependent  solution  for  simple  diffusion  with  a  constant  source. 


in  program  REDSTONE  and  in  the  analytic  solution  results  in  the  comparison 
shown  in  fig.  3.  Exact  agreement  is  obtained  by  REDSTONE  as  would  be  expected 
for  a  strictly  linear  smoothly  varying  solution. 
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8.1.2  Time- Dependent  Solution 

The  problem  that  was  solved  above  was  simple.  A  more  difficult  problem 
involves  the  time  dependence  since  this  will  introduce  near  discontinuities  in  the 
slope  of  the  solution  at  the  grain  boundary.  The  time  dependent  equation  for  a 
single  bubble  group  in  1-D  spherical  geometry  is 


dC  D  d  (  n  dC\  „ 
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The  same  boundary  conditions  are  used  as  before  with  the  added  initial  condition  of 
C  =  0  at  t=0.  For  these  conditions,  a  series  solution  can  be  obtained  [6],  i.e. 
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Fig.  3  Linear  Steady-State  Comparison  with  Analytic  Solution 


119 


K  (  2  2\  ^  njEr\  ( 


Dn^Tp't\ 

K  ] 


(  203  ) 


Calculations  using  REDSTONE  for  time  =  3x10^  secs  are  compared  with  the  series 
solution  in  fig.  4.  Although  the  REDSTONE  calculation  represents  the  series 
solution  well  near  the  center  of  the  grain,  it  oscillates  wildly  at  the  outside  of  the 
grain.  This  oscillation  is  caused  by  the  inability  of  the  interpolation  fimctions  to 
adequately  represent  the  behavior  of  the  concentration  in  a  region  of  sharp  change 
and  is  similar  to  shock  behavior  in  mechanical  modeling.  Ideally,  one  would  want 
interpolation  functions  which  followed  the  profile  of  the  concentration  exactly.  This 
is  not  generally  possible  for  real  world  problems.  However,  one  can  improve  the 
approximation  of  the  concentration  profile  in  several  ways,  two  which  are  used  here. 

The  first  improvement  involves  adjusting  the  interpolation  functions  in  the 
outermost  element  to  better  approximate  the  sharp  decline  in  concentration  at  the 
grain  boundary.  Burnett  [120]  notes  that  the  slope  of  three  node  interpolation 
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Fig.  4  Linear  Transient  Comparison  with  Series  Solution 
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functions  at  one  edge  of  an  element  can  be  made  infinite  by  moving  the  central 
element  off  center  to  a  position  0.75  of  the  distance  between  the  two  outer  nodes. 
This  arrangement  and  the  effect  on  the  interpolation  functions  are  shown  in  fig.  5. 
The  effect  on  the  results  from  REDSTONE  are  compared  with  the  original 
REDSTONE  calculation  in  fig.  6.  A  noticeable  lessening  of  oscillatory  behavior  has 
occurred.  However,  oscillations  remain  because  the  interpolation  functions  still  do 
not  exactly  represent  the  concentration. 

A  second  way  to  improve  the  representation  is  to  increase  the  number  of 
elements  and  nodes  in  the  vicinity  of  the  change.  This  was  tried  for  REDSTONE 
using  the  element  scheme  as  outlined  in  fig.  7.  Although  the  same  number  of 
elements  as  used  previously  are  present,  their  spacings  have  changed  so  that  much 


Fig.  5  Change  in  Interpolation  Functions  with  Middle  Node  Shift 
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Fig  6  Comparison  of  Results  with  and  without  Node  Shift 


Even  Spacing  Graduated  Spacing 


Fig.  7  Graduated  Node  Density  Spacing  in  Fuel  Grain 
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more  detail  can  be  modeled  near  the  grain  boundary.  The  elements  sizes  gradually 
diminished  from  inside  the  grain  to  the  surface.  The  REDSTONE  calculation  using 
this  elemental  scheme  (and  the  node  shift)  is  compared  to  the  previous  REDSTONE 
calculations  in  fig.  8.  Oscillatory  behavior  has  now  become  nearly  negligible,  being 
reduced  by  a  factor  of  at  least  1 000.  Increasing  the  nodal  density  in  the  region  of 
change  has  caused  a  significant  improvement  in  REDSTONE’s  ability  to  model 
actual  distribution  profiles. 

It  is  apparent  from  the  above  analysis  that  a  combination  of  correct  finite 
element  programing  and  an  intelligent  approach  to  the  problem  are  required  to 
achieve  acceptable  results.  Through  the  methods  presented  above,  REDSTONE  has 
been  shown  to  be  capable  of  adequately  representing  the  problem  at  hand. 

8.2  Validation  of  Program  REDSTONE  Results 

A  program  may  work  well  yet  still  provide  unacceptable  results  if  the  correct 
equations  are  not  programmed.  Comparison  with  known  experimental  results 
provides  indications  of  how  well  a  program  models  the  real  world.  Those 
comparisons  for  REDSTONE  are  presented  here. 

8.2. 1  Additional  Constraints 

Initial  calculations  of  multi-bubble  problems  utilized  the  full  radius 
dependent  formulations  of  bubble  radii  change  presented  in  Chapter  VII.  It  was 
quickly  realized  that  the  bubble  radii  calculations  were  extremely  sensitive  to 
changes  in  the  radii.  If  bubble  radii  changed  too  much  in  a  given  time  step,  the 
program  would  begin  to  diverge.  Substitution  of  the  simpler  partial-radius 
dependent  formulation  improved  the  stability  of  the  program  by  slowing  the  rate  of 
change.  However,  the  bubble  radii  calculation  was  still  very  sensitive  to  the  size  of 
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change.  A  contributing  factor  to  program  instability  was  found  in  the  calculation  of 
bubble  pressures.  Due  to  the  form  of  the  reduced  van  der  Waals  equation  of  state,  it 
was  possible  to  calculate  negative  pressures  if  the  radii  became  too  small.  For  the 
small  highly  pressurized  bubbles  here,  this  occasionally  happened  in  the  course  of 
calculations.  These  negative  pressures  affected  calculations  through  the  vacancy 
emission  terms  in  both  vacancy  concentration  and  radii  calculation.  The  bubble 
pressure  was  found  to  be  very  sensitive  to  bubble  radius  in  this  region  and  the 
overly-small  radii  were  always  only  a  few  percent  away  from  a  radius  that  would 
produce  a  positive  pressure.  This  problem  was  solved  by  limiting  the  pressure  of  a 
bubble  to  a  maximum  pressure.  The  yield  pressure  is  a  natural  choice  for  this 
maximum  pressure.  Unfortunately,  the  yield  pressure  of  UN  as  measured  in  the 
laboratory  exhibits  strain  hardening,  making  if  difficult  to  ascertain  what  the 
appropriate  value  to  use  within  the  grain  is.  Because  of  this  a  crude  approximation 
of  the  yield  pressure  was  made  by  using  the  curves  of  Werner  and  Blank  [84],  as 
mentioned  in  Chapter  HI.  The  hope  was  that  these  curves  could  at  least  provide  a 
temperature  dependence  and  an  indication  of  the  pressure  needed  for  plastic 
deformation.  Fortunately,  the  bubble  pressure  is  very  sensitive  to  radius  in  this 
region.  Conversely,  the  radius  is  insensitive  to  pressure.  Thus  large  changes  in 
pressure  do  not  significantly  affect  the  radius  in  the  region  near  the  yield  point.  A 
few  preliminary  calculations  using  REDSTONE  proved  this  to  be  the  case.  Thus  the 
correlation  presented  in  Chapter  III  was  successfully  used  to  eliminate  problems 
with  negative  calculated  pressures. 

This,  however,  did  not  eliminate  the  sensitivity  of  the  radii  calculations  to 
changes  in  bubble  radii,  especially  during  the  severe  transient  to  be  modeled  here.  It 
was  found  that  the  vacancy  emission  term  was  still  the  destabilizing  term.  The 
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unique  formulation  of  this  term  due  to  the  pressure  dependence  may  have  caused  a 
loss  of  too  much  of  the  non-linearities.  This  could  have  then  resulted  in  the  inability 
of  the  correction  term  to  drive  the  solution  in  the  correct  direction  when  changes 
became  too  large.  It  was  found  that  the  calculations  could  be  stabilized  by  limiting 
the  maximum  size  of  the  bubbles  to  the  equilibrium  size.  This  was  an  acceptable 
constraint  since  bubbles  would  be  expected  to  remain  at  or  below  their  equilibrium 
size  in  steady-state  or  ramp-up  conditions.  However,  in  a  ramp-down  condition, 
bubbles  are  likely  to  become  underpressurized  resulting  in  sizes  larger  than  the 
equilibrium  size.  Enforcing  a  maximum  bubble  size  constraint  in  this  situation 
causes  the  bubbles  to  attain  their  equilibrium  size  which  might  affect  the  calculated 
results. 

8.2.2  Comparison  with  Experimental  Studies 

With  program  calculations  stabilized,  REDSTONE  was  ready  to  compare 
with  experimental  studies.  Unfortunately,  the  number  of  experimental  studies  using 
UN  are  a  small  group.  In  addition,  those  studies  measuring  the  types  of  data  needed 
for  comparison  to  intragranular  code  predictions  are  even  fewer.  No  studies  of  UN 
were  found  where  a  microscopic  determination  of  swelling  components  were  made. 
Studies  which  were  found  included  only  those  measuring  macroscopic  swelling  and 
gas  release.  Most  of  these  studies,  especially  the  more  recent  ones,  have  dealt  with 
fuel  temperatures  much  higher  than  the  those  of  interest  here.  This  further  reduces 
the  candidates  for  comparison  since  program  REDSTONE  has  the  added  restriction 
that  the  species  modeled  cannot  be  too  mobile.  Large  bubbles  are  produced  when 
the  species  are  mobile  enough  for  significant  coalescence  to  occur,  as  with  high 
temperatures.  This  is  a  restriction  because  only  individual  bubble  sizes  are  modeled 
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rather  than  bubble  groups  with  averaged  properties.  As  a  result,  large  bubbles  can 
only  be  modeled  using  very  large  program  arrays,  which  take  up  large  amounts  of 
computer  storage  and  exponentially  increase  the  running  times.  Local  computer 
constraints  restricted  the  size  of  the  problem  here  to  50  bubble  sizes.  Even  at  50 
bubble  sizes,  the  running  times  increased  over  six-fold  when  the  number  of  bubble 
sizes  was  doubled.  Thus,  only  studies  utilizing  lower  temperatures,  shorter  bumups, 
or  both  were  candidates  for  comparison. 

When  comparing  REDSTONE  results  to  experimental  release  and  swelling 
studies,  one  should  keep  in  mind  that  these  studies  are  based  on  release  from  fuel 
pellets,  not  individual  fuel  grains.  In  general,  since  gas  released  from  the  grain  will 
be  collected  on  grain  boundaries  before  release  from  the  pellet,  the  release  reported 
in  the  studies  should  be  lower  than  that  predicted  by  an  intragranular  code. 

However,  some  minimum  amount  of  release  would  be  expected  from  secondary 
knockouts  and  recoil  atoms  on  the  pellet  surface  and  in  any  open  porosity.  Direct 
recoil  release  will  come  from  fissions  within  approximately  lOpm  of  the  pellet 
surface,  while  knock-ons,  with  less  energy,  will  only  come  from  about  200 A  deep 
into  the  pellet  [5].  Olander  derives  theoretical  equations  for  release  rates  for  both 
these  mechanisms.  After  noting  that  the  equations  only  predict  a  6%  release  from 
10%  bumup  fuel  in  a  0.6  diameter  pellet,  he  concludes  that  these  mechanisms  are 
not  very  important  in  reactor  fuel-element  performance.  However,  this  calculated 
release  is  equivalent  to  0.6%  release  per  atom%  bumup,  which  although  low,  is 
significant  when  comparisons  are  to  be  made.  While  Olander’s  result  is  only  a 
theoretical  release  rate,  it  does  indicate  that  very  low,  but  measurable,  releases  of 
fission  gas  from  fuel  could  be  explained  entirely  by  recoil  and  knockout 
mechanisms. 
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With  regards  to  swelling,  it  needs  to  be  noted  that  swelling  at  the  grain 
boundaries  and  from  non-gaseous  fission  products  are  not  included  in 
REDSTONE’S  calculations.  Olander  [5]  has  stated  that  this  additional  volume 
swelling  amounts  to  between  0. 15  and  0.45%  per  atom%  bumup  for  UO2,  depending 
on  the  stoichiometry  of  the  fuel.  A  similar  range  can  be  assumed  for  UN. 
Additionally,  any  significant  release  to  the  grain  boundary  will  cause  larger  swelling 
rates  than  in  the  interior  of  the  grain.  Thus,  swelling  predictions  by  REDSTONE 
should  be  less  than  in  experimental  studies,  with  better  agreement  obtained  for  low 
grain  boundary  release  situations. 

Two  studies  were  found  which  satisfied  the  lower  temperature,  lower  bumup 
criteria  One  was  a  swelling  study  performed  at  the  National  Aeronautical  and  Space 
Administration  (NASA)  in  1973  [121].  The  second  was  a  release  and  swelling  study 
performed  at  the  Jtq^an  Atomic  Energy  Research  Institute  (JAERI)  in  1977  [122].  In 
addition.  Storms  [123]  has  correlated  experimental  release  data  to  theoretical 
density,  bumup,  and  temperature  for  releases  between  0. 1  and  1 5%.  He  compared 
his  correlation  with  similar  correlations  by  Ross  [124],  Baars  [125],  and  Thomas 
etal  [126].  He  notes  that  his  correlation,  like  others,  tends  to  underpredict  release 
for  large  releases.  Although  he  does  not  discuss  it,  his  correlation  also  overpredicts 
release  for  low  bumup,  low  temperature  situations  as  is  the  case  here.  This  is  also 
tme  for  the  correlation  by  Thomas  et  al,  but  not  tme  for  those  by  Ross  and  Baars. 
Partly  this  is  because  there  is  a  large  scatter  in  the  data,  and  the  various  sets  of  data 
used  by  each  were  slightly  different.  However,  Storms  also  discarded  from  his  data 
all  releases  below  0.1%  regardless  of  irradiation  temp)erature.  This  eliminated  most 
of  the  low  temperature,  low  bumup  data,  i.e  85%  of  the  data  below  1300°K  /  3% 
bumup,  and  70%  of  the  data  below  1 500°K  /  3%  bumup.  Additionally,  he  includes 
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data  from  U-Pu-N,  as  does  Thomas  et  al.,  which  typically  showed  much  larger 
It. cases  at  the  same  temperatures  and  bumups.  As  a  result,  these  latter  two 
correlations  are  skewed  to  predict  larger  releases  than  experiment  has  shown  in  this 
low  temperature,  low  bumup  region.  It  should  be  noted,  however,  that  the  amount 
of  error  introduced  by  these  correlations  is  not  generally  important  since  the  releases 
are  low.  In  consideration  of  the  above,  rather  than  use  a  correlation  to  compare  with 
REDSTONE  results,  the  data  tabulated  by  Storms  will  be  used  to  indicate  the  range 
of  release  expected  experimentally  for  each  of  the  studies  below. 

8.2.3  Comparison  with  NASA  Results 

The  study  at  NASA  measured  swelling  of  six  UN  fuel  pins  containing  four 
pellets  each  which  were  irradiated  for  ^proximately  4000  hours  to  bumups  of  up  to 
two  percent.  The  average  cladding  surface  tempjerature  was  1 100°K.  However, 
there  was  some  degree  of  thermal  cycling  due  to  reactor  reloads.  Swelling  data  was 
obtained  by  removing  three  pellets  from  two  of  the  pins  for  measurement.  Of  these, 
one  had  significantly  lower  bumup  which  allowed  the  results  to  be  compared  with 
REDSTONE  calculations.  Reported  fuel  characteristics  and  operational  conditions 
for  this  pellet  are  given  in  table  5. 

The  NASA  swelling  results  are  compared  with  the  REDSTONE  result  in 
table  6.  Also  in  table  6  are  REDSTONE  release  data  and  Storms’  data  for  bumups 
between  0.74  and  1.36%  and  temperatures  below  1400°K.  In  performing  the 
REDSTONE  calculation,  the  fuel  grain  size  and  exact  operating  history  were  not 
known.  The  grain  size  used  was  20nm  diameter  and  no  thermal  cycling  was 
considered. 
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Table  5 


Fuel  Pellet  Characteristics  of  NASA  Study 


Outside  Diameter  (cm) 

0.3807 

Length  (cm) 

0.6363 

Enrichment  (%) 

10 

Theoretical  Density  (%) 

94.3 

Power  Density  (W/cm^) 

724 

Average  Temperature  (°K) 

1222+28 

Bumup  (%) 

0.948 

Table  6 


Com£arisonofNASA_R^ults^with^RE]DSTONE^esi^ 


%  AV/V 

%  Release 

NASA 

0.99-1.8 

— 

Storms 

— 

0.01-0.10 

REDSTONE 

0.367 

0.0567 

The  values  obtained  by  REDSTONE  are  lower  than  other  values  in  table  6. 
Some  of  this  difference  might  be  due  to  the  thermal  cycling  in  the  operating  history. 
Conceivably,  the  thermal  cycling  could  have  involved  short  periods  of  higher 
temperature  causing  greater  coalescence  of  bubbles,  greater  release,  or  intra-pellet 
cracking  which  could  have  reduced  the  overall  density  of  the  pellet  and  conversely 
increased  volume.  Additionally,  the  chemistry  of  the  fuel  pellets  was  not  controlled 
to  ensure  that  stoichiometric  fuel  was  used.  It  was  previously  mentioned  that  other 
investigators  were  having  trouble  with  non-stoichiometry  around  this  period,  which 
would  have  resulted  in  increased  diffusional  properties  [123]. 

Some  additional  swelling  was  almost  certainly  contributed  by  non-gaseous 
fission  products.  When  this  additional  swelling  is  considered,  program 
REDSTONE’S  swelling  result  falls  close  to  the  range  experimentally  observed.  The 
residual  difference  can  be  attributed  to  the  factors  mentioned  above. 
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REDSTONE’S  gas  release  result  falls  in  the  middle  of  the  data  tabulated  by 
Storms.  One  would  expect  release  from  the  grain  to  exceed  that  from  the  fuel  pellet, 
but  in  fact  a  small  intragranular  release  may  not  be  released  outside  the  fuel  pellet  at 
all.  It  is  likely  that  the  smallest  release  results  tabulated  by  Storms  involved  release 
by  the  recoil  and  knockout  mechanisms,  not  by  release  from  grains  in  the  interior  of 
the  fuel. 

8.2.4  Comparison  with  JAERI  Results 

The  investigators  at  JAERI  studied  carbide,  carbonitride,  and  nitride  fuels 
irradiated  to  bumups  of  24000  MWD/t  at  temperatures  between  1264°K  and 
1634°K.  Controllers  were  used  to  maintain  temperatures  within  +50°K  during 
irradiation.  Three  UN  pellets  were  measured  for  both  gas  release  and  swelling. 

Two  of  these  had  low  enough  bumups  to  be  compared  with  REDSTONE  results. 
Reported  fuel  characteristics  and  operational  conditions  for  these  pellets  are  given  in 
table  7. 

The  JAERI  swelling  results  are  compared  with  the  REDSTONE  results  in 
table  8.  Also  in  table  8  are  Storms’  data  for  bumups  equivalent  to  that  in  the  JAERI 
study,  between  1.36  and  2.34%  and  temperatures  below  1400°K.  In  performing  the 
REDSTONE  calculation,  the  fuel  grain  size  used  was  20pm  diameter. 

As  with  the  NASA  study,  the  REDSTONE  results  are  on  the  lower  end  of 
the  swelling  range  observed  experimentally,  although  slightly  better  agreement  is 
obtained.  When  non-gaseous  swelling  is  considered  the  REDSTONE  results  fall  in 
the  middle  of  the  observed  swelling  ranges.  However,  experimental  swelling  data 
can  be  seen  to  vary  widely.  This  variation,  besides  making  comparisons  with 
REDSTONE  results  tenuous,  indicates  that  in  low  fission  gas  swelling  situations. 
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Table  7 


Fuel  Pellet  Characteristics  of  JAERI  Study 


Pellet  1 

Pellet  2 

Outside  Diameter  (cm) 

0.325 

0.345 

Length  (cm) 

0.7 

0.7 

Enrichment  (%) 

20 

5 

Theoretical  Density  (%) 

88-95 

88-95 

Power  Density  (W/cm^) 

1998 

2224 

Average  Temperature  (°K) 

1169 

1344 

Bumup  (MWD/t) 

12000 

16000 

Table  8 

Comparison  of  JAERI  Results  with  REDSTONE  Results 

%  AVA^ 

%  Release 

%  AVA^ 

%  Release 

JAERI  (Pellet  1) 

0.5-1. 8 

1.2 

— 

— 

JAERI  (Pellet  2) 

— 

— 

0.7-2.4 

0.7 

Storms 

— 

0.09-1.91 

— 

0.09-1.91 

REDSTONE 

0.4619 

0.0223 

0.7200 

0.0185 

other  processes  are  contributing  as  much  or  more  to  the  overall  swelling  of  the  fuel. 
Some  of  the  possible  processes  were  discussed  previously.  Based  on  the 
comparisons  made  here,  REDSTONE  appears  to  be  satisfactorily  accounting  for  the 
intragranular  fission  gas  component  in  the  overall  swelling.  In  general,  however,  it 
appears  that  REDSTONE  results  will  be  below  experimentally  observed  swelling 
percentages  by  as  much  as  1.5%. 

Fission  gas  release  results  are  also  compared  in  table  8,  with  REDSTONE 
results  being  below  that  of  the  JAERI  study  but  within  the  experimentally  observed 
range  of  Storms.  Because  the  percentage  release  is  low,  the  same  arguments  apply 
here  as  for  the  NASA  study.  Most  experimentally  measured  release  likely  did  not 
come  from  diffusion  out  of  fuel  grains. 
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8.3  Modeling  of  the  Transient  Scenario 

The  transient  scenario  of  interest  here  involves  sustained  low  temperature 
operation  followed  by  a  rapid  surge  in  power  by  at  least  two  orders  of  magnitude 
which  is  sustained  for  several  minutes.  The  parameters  of  this  scenario  were 
discussed  earlier.  Table  9  lists  those  parameters  which  were  used  in  REDSTONE  to 
model  the  scenario. 

Results  of  the  REDSTONE  calculation  are  tabulated  in  table  10.  These 
results  show  that  the  swelling  before  the  transient  is  small  and  that  the  transient  has 
little  effect  on  the  overall  swelling  and  release,  presumably  because  the  time  at  high 
temperature  is  short.  The  amount  of  gas  released  is  insignificant,  especially  when 
release  from  other  non-diffusional  sources  are  considered.  The  post-transient 
swelling  of  1 .0649%  corresponds  to  a  diametral  increase  of  0.354%.  If  one  adds  to 
this  the  possible  1.5%  swelling  due  to  non-gaseous  processes,  the  diametral  increase 
is  still  only  0.848%.  A  check  of  the  data  in  table  1  will  reveal  that  the  fuel  must 
increase  in  diameter  4.6%  before  coming  into  contact  with  the  liner,  and  assuming 
unrestrained  swelling  against  the  cladding,  must  swell  5.9%  before  causing  design 
failure.  Design  failure  was  specified  as  1%  diametral  increase  of  the  fuel  pin. 
Alternately,  swelling  of  the  fuel  pin  because  of  released  gas  is  not  indicated  since 
the  gas  release  from  the  grain  is  insignificant.  It  is  clear  that  the  results  reported 
here  do  not  indicate  fuel  pin  failure  will  occur  due  to  swelling  or  release  of  fission 
gases  from  the  fuel. 

Specific  characteristics  of  the  solution  are  illustrated  in  figs.  9-12.  Fig.  9 
shows  the  swelling  and  gas  release  during  the  transient  and  high  temperature 
periods.  As  seen  in  fig.  9,  the  swelling  initially  increases  with  the  transient,  lagging 
slightly,  then  continues  to  increase.  The  percentage  release  shows  no  effect  from  the 
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Table  9 


Transient  Scenario  Modeling  Parameters 


Fuel  Characteristics 

Fuel 

UN 

Theoretical  density  (%) 

94.3 

Enrichment  (%) 

93 

Grain  diameter  fum) 

20 

Operational  Characteristics 

Initial  power  (W/cm^) 

52.6 

(equivalent  to  10  kW  electric) 

Initial  temperature  (K) 

500 

Initial  temperature  gradient  (K/cm) 

38 

Time  at  initial  power  (sec) 

3x10* 

Transient  time  (sec) 

30 

Final  power  (W/cm^) 

5260 

(equivalent  to  1  MW  electric) 

Final  temperature  (K) 

1500 

Final  temperature  gradient  (K/cm) 

2550 

Time  at  final  power  (sec) 

1000 

Table  10 

Transient  Scenario  REDSTONE  Results 

Pre-transient 

Post-transient 

Volume  Swelling  (%)  1.0627 

Gas  Release  (%)  * 

1.0649 

5.274x10-5 

*  Pre-transient  eas  release  was  below  the  significance  maintained  by  the  computer. 

transient  and  continues  to  decrease  slightly.  Not  shown  is  the  pre-transient  time 
profile,  but  this  would  have  shown  an  initial  percentage  release  of  gas  which  reaches 
a  maximum  and  then  begins  a  steady  decrease  once  the  bubble  concentration  has 
been  built  up  enough  to  prevent  significant  diffiisional  loss  from  the  grain.  The 
bubbles  do  not  contribute  to  any  significant  release,  either  before  the  transient  or 
after. 

An  inspection  of  the  characteristics  shown  in  the  other  figures  indicates  that 
the  swelling  rate  during  the  high  temperature  region,  although  small,  is  greater  than 
th?it  from  simple  deposition  of  gas  atoms  in  the  lattice.  Fig.  10  shows  the 
concentration  profile  in  the  center  of  the  grain  both  immediately  before  the  transient 
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Fig.  9  Swelling  and  Release  Percentages  After  Onset  of  Transient 
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Fig.  10  Pre-  and  Post-transient  Concentration  Profiles 
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and  at  the  end  of  the  high  temperature  period.  The  concentration  profile 
significantly  changes  only  for  single  gas  atoms  and  in  the  large  bubble  tail  of  the 
profile.  The  single  gas  atom  concentration  is  becoming  depleted,  and  the  larger 
bubbles  are  beginning  to  grow  in  number.  The  highest  concentration  of  bubbles  are 
those  containing  four  gas  atoms  both  before  and  after  the  transient.  It  is  clear  that 
there  is  a  shift  occurring  to  larger  bubbles.  The  question  remains,  however,  on 
whether  this  translates  into  an  increase  in  swelling.  A  check  of  the  bubble  radii 
shown  in  fig.  1 1,  both  before  and  after  the  transient,  shows  little  change  occurring, 
with  only  a  slight  increase  for  larger  bubbles.  This  is  supported  by  fig.  12  which 
shows  the  change  in  bubble  pressures.  The  smaller  bubbles  have  actually  increased 
in  pressure  after  the  transient,  as  would  be  expected  if  the  temperature  rose  but  the 
radii  remained  relatively  constant.  The  largest  bubbles,  on  the  contrary,  started  out 
overpressurized  and  actually  decreased  in  pressure  after  the  transient.  This  could 
only  happen  if  the  radii  had  become  larger.  The  implications  of  the  above  analysis 


Fig.  1 1  Pre-  and  Post-transient  Radii  Profiles 
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is  that  swelling  is  increasing  faster  than  the  deposition  of  new  gas  atoms.  While  the 
changes  in  swelling  are  extremely  small  and  not  important  for  the  scenario  modeled 
here,  they  do  indicate  processes  occurring  which  may  become  of  importance  at 
higher  temperatures.  To  investigate  this  possibility,  the  scenario  above  was  modeled 
at  higher  temperatures.  This  investigation  is  discussed  below. 

8.4  Variation  of  Fuel  Temperature 

An  investigation  into  varying  the  fuel  temperature  was  done  by  increasing 
the  pre-  and  post-transient  temperatures  by  1 00“K,  200°K,  and  300“K.  Thus  the  pre¬ 
transient  temperatures  modeled  were  500"K,  600°K,  700°K,  and  800°K.  The 
corresponding  post-transient  temperatures  were  1500°K,  1600°K,  1700°K,  and 
1800°K.  As  shown  in  fig.  13,  the  swelling  increases  steadily  with  temperature  both 
before  and  after  the  transient.  However,  the  amount  of  increase  is  small  and  not 
important  for  the  scenario  modeled. 


Fig.  12  Pre-  and  Post-transient  Bubble  Pressure  Profiles 
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The  associated  gas  release  presented  in  fig.  14  shows  only  post- transient 
release  as  the  pre-transient  release  was  extremely  small  and  essentially  unaffected 
by  the  temperature.  This  lack  of  pre-transient  temperature  effect  is  supported  by  the 
pre-transient  concentration  profiles  in  fig.  15  which  show  virtually  no  change  in 
profiles  with  temperature.  The  post-transient  release  also  changes  little  with 
temperature  except  for  the  results  at  1 800°K.  The  small  changes  in  release  at 
1500°K,  1600°K,  and  1700°K  indicates  that  bubbles  still  are  not  significantly 
contributing  to  release,  and  in  fact  are  keeping  the  single  gas  atoms  from  being 
released.  The  larger  release  at  1 800°K  may  indicate  that  the  bubbles  are  becoming 
more  mobile  at  this  temperature. 

A  check  of  the  post-transient  concentration  profile  in  fig.  16  supports  this 
assertion.  In  fact,  the  concentration  profiles  at  all  temperatures  show  marked 
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Fig.  13  Volume  Swelling  vs.  Pre-  and  Post-transient  Temperatures 
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Fig.  14  Post-transient  Gas  Release  vs.  Temperature 


Fig.  1 5  Pre-transient  Concentration  Profiles  vs.  Temperature 
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differences  between  each  other.  The  concentration  of  the  largest  bubble  size 
increases  significantly  with  temperature  until,  at  1800°K,  it  appears  to  be  on  the 
leading  edge  of  an  upward  trend  and  is  of  the  same  order  of  magnitude  as  the  highest 
concentration  shown.  Additionally,  the  profiles  at  1700°K  and  1800°K  show 
oscillations  at  the  smaller  bubble  sizes.  This  is  possibly  due  to  the  increasingly 
depleted  concentration  of  single  gas  atoms  with  higher  temperature  which  induces  a 
calculational  shock  to  the  solution.  Increasing  the  number  of  time  steps  in  the 
transient  should  reduce  this  oscillatory  behavior. 

Investigating  the  concentration  profiles  at  the  higher  temperatures  further,  it 
appears  that  the  bubbles  are  becoming  too  mobile  for  REDSTONE  to  be  reliable  at 
1800°K  and  possibly  at  ITOO^K.  Fig.  16  would  indicate  that  there  are  much  larger 
bubbles  at  the  same  or  higher  concentrations  at  1800°K.  These  larger  bubbles  have 
the  potential  to  significantly  increase  the  actual  swelling  over  what  REDSTONE  has 
calculated.  Even  when  the  1 800®K  scenario  was  run  with  50  bubble  sizes,  the 


Fig.  16  Post-transient  Concentration  Profiles  vs.  Temperature 
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profile  remained  relatively  flat.  Without  a  reasonable  indication  that  nearly  all  the 
bubbly  sizes  are  accounted  for,  it  is  questionable  to  rely  on  REDSTONE’s  results. 

So  that  some  measurable  indication  could  be  used,  a  figure-of-merit  was 
developed  to  help  discriminate  when  REDSTONE’s  results  might  be  unreliable. 

This  figure-of-merit  was  simply  a  ratio  of  the  concentration  at  the  maximum  in  the 
profile  to  the  concentration  in  the  largest  bubble  size.  This  is  plotted  in  fig.  17.  The 
effect  of  increasing  temperature  on  post-transient  results  can  clearly  be  seen  in  fig. 

1 7.  The  question  now  is  where  to  draw  the  line  below  which  results  are  considered 
to  be  unreliable. 

For  the  purposes  here,  a  liberal  definition  of  a  reliable  REDSTONE  result 
was  chosen.  If  the  figure-of-merit  was  above  100,  the  result  was  considered  reliable. 
Using  this  criteria,  the  results  at  1800°K  would  clearly  be  unreliable,  while  those  at 
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Fig.  17  Pre-  and  Post-transient  Figures-of-Merit 
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17C)0°K,  while  close,  would  be  considered  reliable.  The  previously  mentioned 
modeling  of  50  bubble  sizes  at  1 800°K  resulted  in  an  equally  poor  figure-of-merit  as 
with  30  bubble  sizes. 

The  most  important  observation  of  the  above  c  .  ?sion  is  that  the  swelling 
and  release  behavior  at  1800°K  remains  unknown.  Indications  from  REDSTONE 
are  that  bubble  mobility  becomes  appreciable  at  this  temperature,  leading  to  large 
bubble  sizes  and  possibly  significant  swelling.  This  is  important  since  scoping 
calculations  indicated  a  fuel  centerline  temperature  of  around  1800°K,  even  though 
the  average  temperature  of  the  fuel  would  remain  cool  enough  that  swelling  would 
not  be  a  problem.  However,  the  calculations  here  assume  a  pjerfectly  uniform 
system.  In  an  operating  reactor,  fuel  will  vary  in  temperature  by  position  in  the  core 
and  possible  imperfections  in  the  fuel  or  cladding.  Thus,  some  pellets  might  have 
average  temperatures  near  or  above  1800®K,  and  the  temperatures  would  certainly 
increase  if  final  power  were  boosted  with  minimal  changes  in  design.  Because  this 
temperature  is  close  to  the  temperatures  expected  to  be  present  during  the  current 
transient  scenario,  knowing  the  gas  behavior  at  1800®K  is  an  important  design 
prerequisite.  Unfortunately,  REDSTONE  is  not  readily  enable  of  producing  this 
without  some  sort  of  bubble  grouping  scheme.  Additionally,  in  consideration  of  the 
large  uncertainties  in  some  of  the  input  properties  used  in  REDSTONE,  one  should 
not  discount  the  possibility  that  increased  bubble  mobility  actually  becomes  a 
problem  at  lower  temperatures  than  predicted  by  REDSTONE. 

8.5  Variation  of  Property  Data 

Sensitivity  studies  were  performed  to  determine  whether  the  uncertainties  in 
properties  could  cause  a  significant  effect  on  the  results.  The  propjerties  investigated 
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included  the  xenon,  surface,  and  irradiation  enhanced  diffusion  coefficients,  the 
surface  energy,  and  the  resolution  parameter.  The  same  transient  scenario 
investigated  above  was  used  except  that  a  much  coarser  time  stepping  routine  was 
used  to  speed  program  operation. 

Each  of  the  following  discussions  utilize  a  gas  release  plot  showing  post¬ 
transient  gas  release,  a  pre-  and  post-transient  swelling  plot,  and  the  figure-of-merit 
plot  discussed  previously.  Except  for  surface  energy,  the  properties  were  varied 
through  scaling  by  factors  of  0. 1,  1,  10,  and  100.  The  resolution  parameter  was 
additionally  scaled  by  0.01.  Surface  energy  was  compared  at  780,  1000,  1500,  and 
2000  ergs/cm*. 

8.5. 1  Xenon  Diffusion  Coefficient 

The  gas  release  results  obtained  by  scaling  the  xenon  diffusion  coefficient 
are  presented  in  fig.  18.  As  expected,  an  increase  in  post-transient  release  is  seen 


Scaling  Factor 

Fig.  18  Scaled  Xenon  Diffusion  Coefficient  Gas  Release 
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with  increased  xenon  diffusion  coefficients.  The  reason  for  the  leveling  off  for  the 
highest  diffusivities  is  unknown.  Conceivably,  this  may  have  resulted  from  the 
quick  depletion  of  gas  atoms  through  release  or  capture. 

The  swelling  results  (fig.  19)  show  similar  trends.  As  expected,  the 
pre-transient  swelling  remains  completely  unaffected  by  the  xenon  diffusion 
coefficient  since  irradiation  enhanced  diffusion  is  by  far  the  dominant  process  at  the 
low  pre-transient  temperature  The  post-transient  results  show  the  expected  increase 
in  swelling  as  increased  gas  atom  capture  occurs,  but  with  an  unexplained  dip  in 
swelling  .for  the  highest  diffiisivity.  This  dip  corresponds  with  the  leveling  off  in 
release  results.  A  check  of  the  figure-of-merits  in  fig.  20  indicates  that  all  results 
should  be  reliable. 
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Fig.  19  Scaled  Xenon  Diffusion  Coefficient  Swelling 
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8.5.2  Surface  Diffusion  Coefficient 

The  gas  release  results  obtained  by  scaling  the  surface  diffusion  coefficient 
are  presented  in  fig.  21 .  Since  bubbles  do  not  appear  to  contribute  to  gas  release,  it 
was  expected  that  the  surface  diffusion  would  have  minimal  effect  on  the  gas 
release.  This  is  generally  true  except  for  the  highest  diffiisivity.  At  100  times  the 
surface  diffusion,  it  appears  that  bubbles  began  contributing  to  the  gas  release. 

Again,  however,  the  amount  released  is  extremely  small. 

Because  the  release  is  small,  one  would  expect  that  increased  surface 
diffusion  would  result  in  more  coalescences  leading  to  higher  concentrations  of 
large  bubbles  and  higher  swelling.  Indeed,  this  is  seen  to  be  the  case  in  fig.  22. 
Again,  because  of  the  already  very  small  mobilities  at  the  low  pre-transient 
temperature,  the  pre-transient  bubble  motion  appears  unaffected  by  increased 
surface  diffusion.  A  check  of  the  figure-of-merit  (fig.  23)  supports  this  assertion  and 
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Fig.  20  Scaled  Xenon  Diffusion  Coefficient  Figures-of-Merit 


Scaling  Factor 

Fig.  21  Scaled  Surface  Diffusion  Coefficient  Gas  Release 
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Fig.  22  Scaled  Surface  Diffusion  Coefficient  Swelling 


146 


the  assertion  that  more  coalescences  are  occurring.  However,  the  figure-of-merit  for 
the  post-transient  swelling  at  the  highest  diffiisivity  comes  close  to  the  reliability 
cutoff 

8.5.3  Irradiation  Enhanced  Diffusion  Coefficient 

The  gas  release  results  obtained  by  scaling  the  irradiation  enhanced  diffusion 
coefficient  are  presented  in  fig.  24.  Although  irradiation  enhanced  diffusion  is 
important  in  the  pre-transient  region,  release  results  in  this  region  remained  below 
the  level  of  significance  stored  in  the  computer  for  all  scalings  of  the  diffusion 
coefficient.  When  the  effect  on  swelling  (fig.  25)  is  considered,  it  is  seen  that 
swelling  increases  with  increased  irradiation  enhanced  diffusion  during  the  pre¬ 
transient  period.  This  trend  is  the  result  of  increased  gas  atom  capture  rather  than 
release.  The  post-transient  release  results  suggest  this  same  trend  since  greater  gas 
atom  capture  in  the  pre-transient  period  would  result  in  a  reduced  gas  concentration 
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Fig.  23  Scaled  Surface  Diffusion  Coefficient  Figures-of-Merit 
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being  available  for  release  after  the  transient.  The  bubbles  do  not  appear  to 
contribute  significantly  to  release. 

The  swelling  results  shown  in  fig.  25  indicate  both  pre-  and  post-transient 
swelling  increases  with  increasing  irradiation  enhanced  diffusion.  The  pre-transient 
trend  was  discussed  above.  The  post-transient  swelling  increases  because  of  the 
corresponding  increase  in  the  pre-transient  period.  Interestingly,  however,  the 
post-transient  results  exhibit  smaller  increases  compared  to  the  pre-transient 
increases.  If  bubbles  did  not  have  time  to  move  appreciably  during  the  high 
temperature  phase,  one  might  expect  a  closer  relationship  between  the  two  trends. 
Since  they  appear  different,  this  could  indicate  that  the  post-transient  swelling  is 
only  partially  dependent  on  the  pre-transient  condition,  and  that  coalescences  after 
the  transient  play  a  significant  role  in  setting  up  the  post-transient  profile.  However, 
a  check  of  the  figures-of-merit  (fig.  26)  before  and  after  the  transient  indicate  that 
the  results  at  he  highest  diffiisivity  are  near  the  reliability  limit,  so  there  may  be 
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Fig.  26  Scaled  Enhanced  Diffusion  Coefficient  Figures- of- Merit 
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some  uncertainty  in  the  exact  values  calculated  at  100  times  the  diffusion 
coefficient.. 

8.5.4  Surface  Energy 

The  gas  release  results  obtained  by  setting  the  surface  energy  at  various 
energies  are  presented  in  fig.  27.  Since  the  surface  energy’s  main  effect  is  in  the 
determination  of  the  size  of  the  bubbles,  it  was  not  expected  that  gas  release  from 
the  gram  would  be  affected  by  changing  the  surface  energy.  This  is  demonstrated  in 
the  release  results  obtained. 

The  swelling  results  (fig.  28)  clearly  show  the  effect  of  changing  the  bubble 
sizes.  As  expected,  increasing  the  surface  energy  decreases  the  swelling,  and 
decreasing  the  surface  energy  increases  swelling.  There  does  not  appear  to  be  any 
significant  effect  on  bubble  mobility  since  the  two  curves  in  fig.  28  exhibit  similar 
trends.  A  check  of  the  figures-of-merit  (fig.  29)  indicates  that  all  data  are  reliable. 
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Fig.  27  Surface  Energy  Gas  Release 
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Fig.  28  Surface  Energy  Swelling 
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Fig.  29  Surface  Energy  Figures-of-Merit 
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8.5.5  Resolution  Parameter 

The  gas  release  results  obtained  by  scaling  the  resolution  parameter  are 
presented  in  fig.  30.  A  significant  effect  on  release  can  be  seen.  Since  the 
resolution  parameter  used  here  is  one  to  two  orders  of  magnitude  lower  than 
experimentally  observed  for  UO2,  these  results  show  that  UO2  equivalent  resolution 
results  in  a  four  order  of  magnitude  increase  in  release.  The  amount  of  release  is 
still  not  large  enough  to  be  of  importance.  The  cause  of  the  increase  is  easily 
understood.  As  resolution  increases,  more  gas  atoms  are  kept  free  in  the  lattice. 
These  single  gas  atoms  are  then  more  likely  to  diffuse  out  of  the  grain  than  be 
captured  and  kept  by  a  bubble. 

A  look  at  the  swelling  results  in  fig.  31  show  interesting  trends.  The  pre¬ 
transient  swelling  exhibits  a  maximum  at  a  scaling  factor  of  one  with  decreased 
swelling  as  one  either  increases  or  decreases  the  resolution  rate.  The  decrease  in 
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Fig.  30  Scaled  Resolution  Parameter  Gas  Release 
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swelling  with  increasing  resolution  is  easily  understood.  The  decrease  in  swelling 
with  decreasing  resolution  is  confusing  until  one  checks  the  associated 
figures-of-merit  in  fig.  32.  From  these,  it  is  clear  that  the  concentration  profile  is 
flattening  as  resolution  decreases.  The  lowest  resolution  results  are  clearly 
unreliable,  and  the  next-to-lowest  results  are  suspect.  The  pre-transient  result  at  0. 1 
cannot  be  eliminated  based  on  the  criteria  that  has  been  set  up.  The  reason  for  this 
decrease  is  not  clear.  Since  the  bubbles  are  essentially  frozen  at  the  low 
pre-transient  temperature  and  thus  bubble  concentration  grows  primarily  through 
gas  atom  capture,  one  could  argue  for  some  decrease  in  swelling  with  lower 
resolution  if  the  single  gas  concentration  became  so  low  that  few  were  left  to 
adequately  develop  a  significant  concentration  of  large  bubbles.  However,  this 
argument  seems  tenuous,  at  best. 

Considering,  now,  the  post-transient  swelling  one  sees  that  it  exhibits  a  much 
milder  dependence  on  resolution.  As  expected,  it  decreases  with  increasing 
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Fig.  3 1  Scaled  Resolution  Parameter  Swelling 


resolution  rate.  However,  it  appears  that  this  is  further  evidence  that  post-transient 
mobilities  may  have  had  an  effect  in  setting  up  the  post-transient  concentration 
profile. 


0.001  0.01  0.1  1  10  100 


Scaling  Factor 


Fig.  32  Scaled  Resolution  Parameter  Figures-of-Merit 
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CHAPTER  IX 

SUMMARY  AND  CONCLUSIONS 


9.1  Summary 

In  order  to  investigate  whether  fission  gas  swelling  and  release  would  be 
significant  factors  in  a  space  based  nuclear  reactor  operating  under  the  Strategic 
Defense  Initiative  (SDI)  program,  a  finite  element  program  was  developed  to  model 
the  transient  fission  gas  behavior  in  a  single  uranium  nitride  fuel  grain.  Only  the 
intragranular  problem  was  investigated  because  of  the  importance  of  the 
intragranular  fission  gas  as  a  precursor  to  all  other  fission  gas  based  processes,  and 
because  it  was  necessary  to  limit  the  investigation  to  a  reasonably  sized  problem. 

The  program  developed  was  named  REDSTONE  (Routine  for  Evaluating 
Dynamic  Swelling  in  Transient  Operational  Nuclear  Environments)  and  was  based 
on  formulation  of  the  governing  I-D,  spherical  geometry  diffusion  equations  using 
the  Galerkin  method.  The  equations  were  incrementalized  so  that  non-linear 
dependences  would  be  retained,  and  the  solution  was  obtained  through  iteration. 
Calculations  were  divided  into  three  sets  of  governing  equations:  bubble 
concentration  changes,  bulk  vacancy  and  interstitial  concentration  changes,  and 
bubble  radii  changes.  The  bubble  concentration  and  radii  equations  were  derived 
from  considerations  of  production  and  loss  mechanisms  of  individual  bubble  groups 
defined  by  the  number  of  gas  atoms  in  the  bubble.  This  is  in  contrast  to  similar 
investigations  for  UOj  which  have  utilized  groups  with  average  properties 
describing  bubbles  containing  a  range  of  gas  atoms.  The  choice  to  model  individual 
bubbles  placed  limitations  on  the  types  of  scenarios  REDSTONE  could  model.  The 
limitation  was  due  to  the  increased  storage  requirements  and  running  times  needed 
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to  model  a  large  range  of  bubble  groups.  Since  higher  temperatures  and  longer 
bumups  cause  larger  bubbles  to  form,  REDSTONE’s  calculations  were  limited  to 
those  scenarios  where  low  temperatures,  low  bumups,  or  both  were  present. 

The  bubble  concentration  equations  incorporated  the  production  and  loss 
mechanisms  of  diffusion,  production,  resolution,  homogeneous  nucleation,  random 
gas  atom  capture,  biased  gas  atom  capture  (due  to  thermal  gradients),  random 
coalescence,  and  biased  coalescence.  In  addition,  effects  due  to  biased  diffusion  out 
of  the  grain  were  included  in  a  separate  routine  which  calculated  the  biased  release 
by  the  overlap  of  two  spheres  shifted  by  an  amount  consistent  with  the  bubble 
velocity  and  the  time  period. 

The  bulk  vacancy  and  interstitial  concentration  equations  incorporated  the 
production  and  loss  mechanisms  of  diffusion,  production,  thermal  emission  from 
bubbles  and  dislocations,  capture  by  bubbles  and  dislocations,  and  recombination. 

The  bubble  radii  equations  were  derived  from  considerations  of  the  number 
of  vacancies  and  interstitials  entering  and  leaving  a  bubble.  Governing  processes 
from  the  bubble  concentration  equations  were  incorporated  and  modified  by  the 
associated  vacancies  involved  during  the  interactions.  In  addition,  direct  interaction 
with  vacancies  and  interstitials  was  incorporated  through  appropriate  terms  for 
thermal  emission  of  vacancies  and  capture  of  vacancies  and  interstitials.  Originally, 
every  radius  dependent  term  in  these  equations  was  updated  during  iterations. 
However,  instabilities  were  encountered  which  were  lessened  when  the  terms 
coming  directly  from  the  concentration  equation  were  held  constant  during 
iterations.  Holding  these  terms  constant  did  not  eliminate  all  the  causes  of 
instability,  and  further  restrictions  had  to  be  imposed  on  bubble  radii  calculations. 
The  first  of  these  was  the  enforcement  of  a  minimum  radius  a  bubble  could  have. 
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This  was  necessary  because  of  the  occasional  calculation  of  negative  pressures 
based  on  the  radius.  The  minimum  radius  was  determined  based  on  the  yield 
pressure  of  uranium  nitride.  The  yield  pressure  was  determined  from  a  crude  least 
squares  fit  to  available  stress-strain  curves.  A  second  restriction  imposed  was  the 
enforcement  of  a  maximum  radius  the  bubble  could  have.  Tfre  maximum  radius  was 
the  equilibrium  radius  of  the  bubble.  This  avoided  a  stability  problem  in  the  thermal 
emission  term,  presumably  caused  by  a  loss  of  too  much  non-linearity  in 
formulating  the  correction  part  of  this  term.  These  two  additional  restrictions 
stabilized  the  radii  calculations,  but  added  additional  constraints  on  what 
REDSTONE  could  model.  Steady-state  or  ramp-up  scenarios  would  still  be  able  to 
have  their  overpressurized  non-equilibrium  bubbles  modeled.  However,  ramp-down 
scenarios  would  have  their  underpressurized  bubbles  modeled  with  equilibrium 
bubbles. 

The  validity  of  REDSTONE  calculations  were  checked  against  analytical 
solutions  for  internal  consistency  and  against  experimental  studies  for  agreement 
with  the  swelling  and  release  results.  Internal  consistency  checks  validated  that 
REDSTONE’S  basic  structure  was  correctly  programmed.  These  checks  also 
indicated  that  the  best  agreement  with  transient  results  would  be  obtained  by  using 
an  element  scheme  with  smaller  and  smaller  elements  as  the  grain  surface  is 
approached,  and  by  shifting  the  middle  node  in  the  outermost  element  towards  the 
grain  surface.  The  checks  performed  against  experimental  studies  showed 
REDSTONE  was  satisfactorily  accounting  for  the  intragranular  fission  gas  portion 
of  swelling,  but  that  other  non-gaseous  sources  of  swelling  accounted  for  up  to  1.5% 
additional  swelling. 
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Once  validated,  REDSTONE  was  used  to  calculate  the  swelling  and  gas 
release  associated  with  a  hypothetical  SDI  scenario.  Since  an  SDI  nuclear  system 
has  not  been  specified,  the  reactor  parameters  used  here  were  obtained  from  the 
SP- 1 00  nuclear  reactor.  This  reactor  has  a  design  failure  point  of  1  %  diametral 
increase  of  the  fuel  pin.  Operating  conditions  for  the  scenario  were  calculated  by  a 
SP- 1 00  systems  code  developed  at  Texas  A&M  University.  The  SDI  scenario 
involved  10  years  of  low  power  operation  at  10  kW  followed  by  a  30  second  rise  to 
1  MW.  The  1  MW  power  was  maintained  for  approximately  15  minutes.  Fuel 
temperature  was  500°K  at  low  power  and  1500°K  at  high  power.  The  fission  gas 
swelling  and  release  results  calculated  by  REDSTONE  indicated  that  swelling 
would  only  be  slightly  above  1%  and  release  would  be  negligible.  Even  when 
swelling  due  to  non-gaseous  compionents  is  added  to  this,  the  diametral  increase  in 
the  fuel  pellet  remains  below  1%,  with  almost  6%  needed  for  design  failure. 

This  same  scenario  was  run  at  higher  temperatures  to  see  if  uncertainties  in 
the  operating  parameters  could  cause  significantly  different  results.  Both  the  low 
power  and  high  power  temperatures  were  increased  by  100°K,  200°K,  and  300°K. 
The  increased  temperatures  had  no  significant  effect  during  the  low  jxjwer  period. 
After  the  high  power  period,  the  swelling  and  release  results  changed  slightly,  but 
the  differences  were  unimportant.  However,  further  investigation  into  the  calculated 
solutions  revealed  that  the  results  at  the  highest  temperature  of  1800°K  appeared 
unreliable  because  high  bubble  mobilities  were  causing  a  shift  to  bubbles  larger  than 
those  modeled  by  REDSTOI'IE.  Thus  the  swelling  calculated  by  REDSTONE  was 
smaller  than  if  all  possible  bubble  groups  had  been  considered.  Whether  the  actual 
swelling  would  have  been  important  for  failure  considerations  is  not  known  since 
REDSTONE  is  not  capable  of  modeling  the  bubble  sizes  required.  REDSTONE’S 
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results  are  important  indications,  however,  since  centerline  fuel  temperatures  were 
expected  to  be  near  1 800°K. 

It  was  apparent  from  REDSTONE’S  results  that  further  investigations  into 
the  swelling  behavior  near  1800°K  were  needed.  Additional  uncertainties  in 
property  values  used  in  REDSTONE  and  in  the  operating  parameters  of  the  fuel 
might  even  lower  this  transition  temperature  to  1700°K  or  lower. 

Many  of  the  property  values  used  in  REDSTONE  had  to  be  gleaned  from  a 
small  library  of  studies.  In  some  cases,  only  one  value  may  have  been  available  for 
use.  In  other  cases,  the  values  had  to  be  taken  from  other  fuels  such  as  UC  or  UOj. 
As  a  result,  most  of  the  properties  have  a  great  deal  of  uncertainty  associated  with 
them.  In  an  effort  to  understand  the  effect  these  uncertainties  might  have  on  the 
solution,  several  of  the  properties  were  varied  and  the  results  compared.  These 
included  the  xenon,  surface,  and  irradiation  enhanced  diffusion  coefficients,  the 
surface  energy,  and  the  resolution  parameter.  Except  for  surface  energy,  the 
properties  were  varied  by  two  orders  of  magnitude  on  either  side  of  the  programmed 
value.  Most  of  the  changes  caused  by  varying  the  properties  were  small  and 
unimportant  with  regards  to  failure.  However,  each  property  showed  definite  trends 
which  might  become  important  at  higher  temperatures.  Thus,  knowing  these  trends 
is  an  important  design  consideration. 

The  largest  effects  were  observed  from  changes  in  irradiation  enhanced 
diffusion  and  resolution.  Although  the  changes  caused  by  varying  these  properties 
were  small,  both  had  solutions  which  exhibited  significant  shifts  to  larger  bubble 
sizes.  The  irradiation  enhanced  diffusion  coefficient  produced  such  results  above 
100  times  the  programmed  value.  The  resolution  parameter  produced  such  results 
below  0.01  of  the  programmed  value.  However,  the  resolution  parameter  used  in 
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REDSTONE  is  already  lower  than  that  experimentally  observed  for  UOj,  so  it  is 
unlikely  that  the  value  used  in  REDSTONE  is  two  orders  of  magnitude  too  high.  On 
the  contrary,  the  uncertainties  in  the  irradiation  enhanced  diffusion  coefficient  would 
be  more  likely  to  range  towards  larger  values  in  the  same  direction  as  the  fall  in 
reliability. 

In  situations  where  the  bubble  sizes  remain  small,  the  largest  effect  of 
uncertainties  in  property  values  could  come  from  the  resolution  parameter.  The 
value  used  in  REDSTONE  is  one  to  two  orders  of  magnitude  lower  than  in  UOj 
because  it  was  assumed  that  UN  resolution  would  occur  via  a  single  knock-on 
mechanism  versus  a  thermal  spike  mechanism  in  UO.,.  If,  in  fact,  resolution  is  the 
same  in  both  fuels,  then  the  parameter  study  here  indicates  a  reduction  in  the 
swelling  and  a  corresponding  increase  in  release. 

Both  the  irradiation  enhanced  diffusion  and  resolution  property  studies 
suggested  that  the  bubble  profile  after  the  power  rise  was  determined  by  a 
combination  of  the  pre-rise  profile  and  the  mobility  of  the  bubbles  during  the  high 
temperature  period.  Even  though  the  time  at  high  temperatures  was  short,  it 
appeared  that  the  conditions  after  the  power  rise  were  less  a  function  of  the  pre-rise 
conditions  than  hypothesized. 

9.2  Need  for  Further  Work 

Throughout  this  study,  various  areas,  both  external  and  internal  to  this  study, 
were  found  which  deserve  further  investigation.  External  to  the  study,  it  was  clear 
that  information  on  microscopic  UN  properties  is  severely  lacking  and  much  of  it 
dates  to  the  early  1970’s  and  1960’s.  Most  of  the  information  available  is  scattered 
in  the  literature  and  not  very  well  corroborated.  This  includes  data  on  diffusion  of 
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gaseous  species,  the  number  of  vacancies  associated  with  gaseous  species,  the  effect 
of  stoichiometery,  the  effect  of  fission  damage,  surface  energies,  and  nuclear 
processes  such  as  resolution  mechanisms,  fission  damage  production  rate,  and 
nucleation  mechanisms.  These  are  clear  candidates  for  additional  work.  In  addition, 
microscopic  studies  of  swelling  in  UN  needs  to  be  performed  so  that  an 
intragranular  code  can  be  benchmarked  directly  against  the  data  it  produces. 

Internal  to  this  study,  additional  work  was  identified  both  by  additional 
questions  which  arose  and  by  deficiencies  which  were  inevitably  incorporated.  The 
additional  questions  were  ‘what  if  questions  dealing  with  the  scenario  modeled. 

The  study  here  looked  at  a  single  hypothetical  scenario  and  a  few  of  the  parameters 
which  could  affect  it.  There  are  of  course  other  possible  scenarios,  all  which  could 
be  investigated.  These  other  scenarios  might  involve  different  power  histories  as 
well  as  other  reactor  types.  Some  of  the  questions  cannot  be  answered  until  firmer 
design  specifications  for  an  SDI  system  are  developed.  Additionally,  sensitivity 
studies  such  as  the  ones  performed  on  property  uncertainties  remain  areas  for 
additional  work.  Besides  additional  property  data,  calculational  factors  such  as  time 
step  size  and  element  spacing  could  be  investigated.  These  latter  studies  are  useful 
for  fully  characterizing  the  capabilities  of  the  program. 

The  deficiencies  needing  additional  work  deal  with  limitations  on  the  type  of 
scenarios  which  could  be  modeled.  The  main  limitation  was  that  REDSTONE 
needed  to  model  each  individual  bubble  size.  This  prevents  the  calculation  of  results 
in  high  swelling  situations.  As  a  result,  first  on  the  list  of  additional  work  needed  on 
REDSTONE  is  the  incorporation  of  a  bubble  grouping  scheme  which  can  model  the 
larger  bubbles.  The  next  most  important  area  of  improvement  in  REDSTONE  is  the 
calculation  of  bubble  radii  changes.  As  discussed  previously,  constraints  built  into 
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REDSTONE  limit  is  usefulness  in  modeling  underpressurized  bubbles.  These 
constraints  need  to  be  removed.  In  order  to  accomplish  this  the  formulation  of  the 
bubble  radii  equations  needs  to  be  reworked  to  eliminate  instabilities.  Part  of  this 
would  have  to  involve  a  restatement  of  the  equation  of  state  for  xenon  so  that 
negative  pressures  could  not  be  calculated.  Possible  approaches  to  rework  the  radii 
equations  could  be  as  simple  as  replacing  the  full  equations  with  a  relationship 
which  does  not  incorporate  the  terms  of  the  concentration  change  equation  and  does 
not  require  iteration,  or  as  detailed  as  finding  a  way  to  incorporate  all  the  types  of 
equations  into  one  calculation,  and  keeping  the  non-linearities.  Many  investigators 
have  taken  the  former  approach,  and  this  is  probably  the  easiest  to  accomplish. 

9.3  Conclusions 

Several  conclusions  can  be  drawn  from  this  study.  These  are; 

1)  Non-linear  finite  element  techniques  can  be  used  successfully  to  model 
the  transient  and  non-equilibrium  equations  describing  fission  gas 
behavior  in  nuclear  fuel.  The  most  difficult  part  of  the  finite  element 
formulations  as  performed  here  was  in  dealing  with  the  non-equilibrium 
bubble  behavior. 

2)  There  is  a  lack  of  adequate  experimental  studies  on  microscopic  swelling 
of  UN  to  adequately  benchmark  an  intragranular  fission  gas  routine. 

3)  Within  the  benchmarking  limitations,  the  program  developed  for  this 
study  appears  to  satisfactorily  account  for  the  fission  gas  component  of 
swelling  in  the  grain.  Non-gaseous  component  swelling  of  fuel  could 
add  up  to  1.5%  additional  swelling  at  low  bumups. 
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4)  Many  of  the  property  values  needed  to  model  fission  gas  behavior  in  UN 
either  are  not  known  and  must  be  adopted  from  other  materials  or  have 
large  uncertainties  because  of  the  lack  of  corroborating  studies. 
Uncertainties  in  the  irradiation  enhanced  diffusion  coefficient  and  the 
resolution  parameter  appear  to  have  the  greatest  potential  for  causing 
significant  swelling  or  release  changes  at  higher  bumups  and/or  higher 
temperatures. 

5)  The  SDI  scenario  modeled  here  is  not  predicted  to  produce  significant 
fuel  swelling  or  release  provided  the  temperature  remains  below  1700°K. 
Between  1700°K  and  1800°K,  bubble  mobilities  increase  significantly  so 
that  accelerated  swelling  is  possible.  Uncertainties  present  in  property 
values  and  operating  parameters  could  lower  this  temperature  even  more. 
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INCLUDE  (DIFFCOEF)  •  gas  bubbles  using  spherical  geometry. 
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INCLUDE  (ELKF)  * - - - 
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INCLUDE  (B0LT2MAN)  +  SHPFCT(3)*R(NODE3) 

INCLUDE  (SURFENGY) 


- - - - +  lOOFl  =.  ICONDF(IELEM,N) 

Value  of  R  in  real  coordinotes  of  integration  point  I  ID0F2  =  ICONDF(IELEM,(NBBLSZ+N)) 
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else  alpha  =  ALPHA-OMEGA*DIFFV*CONCIPCN)*CONVIP 

RESOLCMD  -  RESOLP-QFLOATCN)  +  *RADIP(N)/CNR0IP(N)*»2 

ENOIF  BETA  =  BETA  -  2.Q0*OMEGA*DIFFV*CONCIPCN)*CONVIP 

SRCFRC(N)  =  8.40128Q0*PI*OMEGA*RAOIP(N)**2/LTPARM**Z  +  *RADIP(NyCNRDIP(N)**3 

TOTSRC  -  TOTSRC  +  SRCFRCCN)*CONCIP(N) 


ALPHA  -  ALPHA+0MEGA*DlFFV*C0NCIP(N)*C0NCVn(l,2)  *- . - . . 

*RADIPCN)/CNRDIP(N)**Z  *  Zero  constants 

•QEXP(0MEGA/(BLTZMN*TEMP(1))*  * - 

C2.Q0*SRFENG/RADIP(N)-BPRSIP(N)+PRESS(1))) 
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Dim:  i  -  max  number  of  bubble  groups  I  •  Dim:  i  =  max  number  of  materials 

- 4.  *  TMPGR(i)  -  temperature  gradient  for  i'th  material 

REAL*16  RAOMAX(50),  RADMIN(50),  ATOM(50),  VELCUMCS0), VOL 0*4(50)  *  Dim:  i  =  max  number  of  materials 

- 4  ♦  PRESS(i)  -  hydrostatic  pressure  for  i’th  material 
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TELK(i,j)  -  temporary  elemental  stiffness  matrix  I  REAL*16  RNAT(6,6),WEIGHT(6,6) 
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APPENDIX  2 
SAMPLE  INPUT 
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* - + 

*  OUPUT  DECK  ?  -  0:  NO  1:  YES  I 

*  - + 

000  0 

* - + 

♦  HEADER  CARDS,  MAX  80  CHARACTERS  I 

♦  - + 

100  UN  Fuel  Grain  Modelling  Using  Nonlinear  Finite  Element 

Analysis 

101  By  Capt.  Daniel  L.  DeForest,  USAF 

102 

103  Sample  Input  Deck 

106  CONDITIONS 

108  Power  Density  =  700.00  lN/cm3 

109  Temp  =  1200  K 

110  Temp  Grad  =  1000  K/cm 

111  External  Pressure  =  1033  Pa 

112  Time  Length  =  1.0E5  s 

• - + 

*  TYPE  OF  UNITS  1 

•  - - - + 

201  2 

• . . . - - - + 

•  NUMBER  OF  GUASS  INTEGRATION  POINTS  1 

•  . - - - + 

202  3 

*- . - . - - - + 

*  LIST  INPUT  DECK  FLAG  I 

*  - - - - + 

203  0 

— - - - - + 

♦  RESTART  FLAG  I 

*  - + 

204  0 

* - + 

*  NE,NBBLSZ,NNELEM,NSPC,NPRD,MAXITER  I 

♦  - + 

300  10  10  3  10  2  50 

♦ - - - - - ^ 

*  R  COORDINATES  IN  ORDER  OF  NODE  I 

♦  - + 

501  0.0Q0 

502  1.25Q-4 

503  2.5Q-4 

504  3.5Q-4 

505  4.5Q-4 

506  5.5Q-4 

507  6.5Q-4 

508  7.0Q-4 

509  7.5Q-4 

510  8.0Q-4 

511  8.5Q-4 
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512  8.9Q-4 

513  9.3Q-4 

514  9.5Q-4 

515  9.7Q-4 

516  9.8Q-4 

517  9.9Q-4 

518  9.925Q-4 

519  9.95Q-4 

520  9.9875Q-4 

521  10.0Q-4 

* - + 

*  SPECIFIED  CONCENTRATIONS  -  NODE  #,  BBLGP,  VALUE  1 


601  21 

1 

0.0Q0 

602  21 

2 

0.0Q0 

603  21 

3 

0.0Q0 

604  21 

4 

0.0Q0 

605  21 

5 

0.0Q0 

606  21 

6 

0.0Q0 

607  21 

7 

0.0Q0 

608  21 

8 

0.0Q0 

609  21 

9 

0.0Q0 

610  21 

10 

0.0Q0 

* 

— + 

•  GRAIN 

RADIUS 

1 

— + 

701  10.Q-4 

• . - . - . . + 


*  INITIAL  POROSITY,  INITIAL  ENRICHMENT,  CONVERGENCE  LIMH  I 

♦  - + 

901  0.057  0.10  l.Q-4 

. - . + 

*  PERIOD  LENGTH,  #TIME  STEPS,  PRINTOUT  INTERVAL  I 

. . . - - - + 

1001  1.0Q5  5  5 

•  - + 

*  MATERIALS  -  POWER, TEMP, TEMP  GRADIENT, PRESSURE  I 

♦  - + 

1102  700.0  1200.0  1000.0  1033.23  700.0  1200.0  1000.0  1033.23 
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APPENDIX  3 
SAMPLE  OUTPUT 
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1 


RRRRRRR 

EEEEEEEE 

DDDDDDD 

ssssss 

iirmn 

000000 

NN 

NN 

EEEEEEEE 

RRRRRRRR 

EEEEEEEE 

DDDDDDDD 

ssssssss 

TTTTTTTT 

00000000 

NNN 

NN 

EEEEEEEE 

RR 

RR 

EE 

DO 

00 

ss 

TT 

00 

00 

NNNN 

NN 

EE 

RR 

RR 

EEEEEE 

OD 

DO 

sssssss 

TT 

00 

00 

NN  NN 

NN 

EEEEEE 

RRRRRRRR 

EEEEEE 

DD 

DD 

sssssss 

TT 

00 

00 

NN  NNNN 

EEEEEE 

RRRRRR 

EE 

DD 

DD 

ss 

TT 

00 

00 

NN 

NNN 

EE 

RR 

RR 

EEEEEEEE 

DDDDDDDD 

ssssssss 

TT 

00000000 

NN 

NN 

EEEEEEEE 

RR 

RR 

EEEEEEEE 

DDDDDDD 

ssssss 

TT 

000000 

NN 

NN 

EEEEEEEE 

— =  PROGRAM  REDSTONE  t  ..n,  ,  .. 

UN  Fuel  Grain  Modelling  Using  Nonlinear  Finite  Element  Analysis 
By  Capt.  Daniel  L.  DeForest,  USAF 


Sample  Input  Deck 
CONDITIONS 

Ponyer  Density  -  700.00  W/cmB 
Temp  -  1200  K 
Temp  Grad  -  1000  K/on 
Externol  Pressure  »  1033  Pa 
Time  Length  -  1.0E5  s 


1 

Total  number  of  elements  -  10 
Total  number  of  nodes  »  21 
Maximum  bubble  groups  per  node  -  10 
Total  number  of  material  types  -  1 
Total  number  of  specified  displacements  =  10 
Total  number  of  time  periods  -  2 
Maximum  number  of  iterations  -  50 
Units  are  in  BRITISH 


Number  of  Guass  integration  points  is  3 


Grain  Diameter 


1.000E-3  microns 
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NODAL  POSITIONS  FROM  CENTER 

RCD 

Global  Nodol  Cloches) 

e.eooeeeiTe 

1.2S0OO0E-4 
2.500e0eE-4 
3.5e00OOE-4 
4.500000E-4 
5.50Oe00E-4 
6.5e0oe0E-4 

7AAAaaac  a 
.  W6WWt"‘H 

7.5O00eeE-4 
8.0000OeE-4 
8.5aOOO0E-4 
8.90e00eE-4 
9.30e000E-4 
9.5e00e0E-4 
9.700000E-4 
9.8e00eoE-4 
9.90O000E-4 
9.925000E-4 
9.95eoe0E-4 
9.987500E-4 
l.e0O0O0E-3 

ELEMENT  CONNECTIVITY 

I - Global  Node  Number - 1 

I - for  the  Local  Node - 1 

I - Numbers  Below - 1 

Element  1234 

T“  r~  3~  T” 

3  4  5  0 

5  6  7  0 

7  8  9  0 


9 

10 

11 

0 

11 

12 

13 

0 

13 

14 

15 

0 

15 

16 

17 

0 

17 

18 

19 

0 

19 

20 

21 

0 

ELEMENT  MATERIAL  PROPERTIES 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 


Material  Initial 
Element>I  Type  ID  Porosity 


5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 

5.7000E-2 
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1  TIME  HISTORY  CONDITIONS 


Time 

Incr 

+ 

Initial 

Power 

Initial 

Temp 

Initial 

Press 

Final 

Power 

Final 

Temp 

Final 

Press 

1 

700. 

1.200E+3 

1.033E+3 

700. 

1.200E+3 

1.033E+3 

2 

0.000E+0 

0.000E40 

0.000E+0 

0.000E-f0 

0.000E40 

0.000E40 

1 

TIME 

STEP  INFORMATION 

Time  Length  Number  of  Time  Step  Printout 

Incr  Cseconds)  Time  Steps  (seconds)  Interval 

+  _  _  _  _  _ 

1  1.000E+05  5  2.00OE-f«4  5 

2  0.000E^  0  -1.000E-f05  0 

1 

SPECIFIED  CONCENTRATIONS 

Global  Node  it  Bubble  Gp  Value 

+  _  _  _ 

21  1  0.000E-f00 

21  2  0.000E400 

21  3  0.000E-f00 

21  4  0.000E400 

21  5  0.000E-^ 

21  6  0.000E-^ 

21  7  0.000E400 

21  8  0.000E-f00 

21  9  0.000E-f00 

21  10  0.000E-f00 

1 

TIME  PERIOD  1 

Time  Step  5 

Time  Step  Size  »  2.0000E+4  seconds 

Time  in  Period  -  1.0000E+5  seconds 

New  Time  -  1.0000E+5  seconds 

Time  step  represents  100*  of  full  time  period. 

Percent  fission  gas  released  -  0.4815  * 

Percent  swelling  -  4.2423E-03% 


Burnup  rate  =  1.129E-02  atom  percent/day 

Burnup  for  period  =  2.614E-03  atom  percent 
Burnup  (total)  -  1.307E-02  atom  percent 
-  122.80  MWD/MT 


Power  density  »  700.0  W/cm3 
Pressure  =  1033.  Pa 

Temperature  =  1200.  K 


228 


NODAL  CONCENTRATIONS 

Bubble  Concentration  Bubble  Radius  Bubble  Pressure 

Croup  Ccni-3)  Ccm)  CPo) 

- Node  ■=  1 - Radius  -  0.000E+00 - 

1  1.907419E+17  2.426457E-08  0.000000E+00 

2  1.152346E+17  3.471774E-08  5.760742E+10 

3  2.970138E+16  3.980550E-08  5.024431E+10 

4  6.723207E+15  4.386735E-08  4.559199E+10 

5  1.396060E+15  4.730536E-08  4.227850E+10 

6  2.723131E+14  5.031653E-08  3.974836E+10 

7  5.061627E+13  5.301388E-08  3.772596E+10 

8  9.051949E+12  5.546897E-08  3.605619E+10 

9  1.577364E+12  5.773034E-08  3.464382E+10 

10  2.973562E+11  5.983260E-08  3.342659E+10 

- Node  -  2 - Radius  -  1.250E-04 - 

1  1.907419E+17  2.426457E-08  0.000000E+00 

2  1.152346E+17  3.471774E-08  5.760742E+10 

3  2.970138E+16  3.980S50E-08  5.024431E+10 

4  6.723206E+15  4.386735E-08  4.559199E+10 

5  1.396059E+15  4.730536E-08  4.227850E+10 

6  2.723131E+14  5.031653E-08  3.974836E+10 

7  5.061627E+13  5.301388E-08  3.772596E+10 

8  9.051950E+12  5.546897E-08  3.605619E+10 

9  1.577365E+12  5.773034E-08  3.464382E+10 

10  2.973562E+11  5.983260E-08  3.342659E+10 

. Node  •  3 - Radius  -  2.500E-04 . 

1  1.907419E+17  2.426457E-08  0.000000E+00 

2  1.152346E+17  3.471774E-08  5.760742E+10 

3  2.970138E+16  3.980550E-08  5.024431E+10 

4  6.723206E+15  4.386735E-08  4.559199E+10 

5  1.396059E+15  4.730536E-08  4.227850E+10 

6  2.723131E+14  5.031653E-08  3.974836E+10 

7  5.061627E+13  5.301388E-08  3.772596E+10 

8  9.051950E+12  5.546897E-08  3.605619E+10 

9  1.577364E+12  5.773034E-08  3.464382E+10 

10  2.973562E+11  5.983260E-08  3.342659E+10 

- Node  «  4 - Radius  ■>  3.500E-04 - 

1  1.907419E+17  2.426457E-08  0.000000E400 

2  1.152346E+17  3,471774E-08  5.760742E+10 

3  2.970138E+16  3.980550E-08  5.024431E+10 

4  6.723206E+15  4.386735E-08  4.559199E+10 

5  1.396059E+15  4.730536E-08  4.227850E+10 

6  2.723131E+14  5.031653E-08  3.974836E+10 

7  5.061627E+13  5.301388E-08  3.772596E+10 

8  9.051950E+12  5.546897E-08  3.605619E+10 

9  1.577365E+12  5.773034E-08  3.464382E+10 

10  2.973562E+11  5.983260E-08  3.342659E+10 

-  Node  «  5  -  Radius  «  4.500E-04 - 


1  1.907419E+17  2.426457E-08  0.000000E+00 

2  1.152346E+17  3.471774E-08  5.760742E+10 

3  2.970138E+16  3.980550E-08  5.024431E+10 
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4 

6.723207E+15 

4.386735E-08 

4.559199E+10 

5 

1.396060E+15 

4.730536E-08 

4.227850E+10 

6 

2,723131E+14 

5.031653E-08 

3.974836E+10 

7 

5.061627E+13 

5.301388E-08 

3.772596E+10 

8 

9.051949E+12 

5.546897E-08 

3.605619E+10 

9 

1.577364E+12 

5. 773034E-08 

3.464382E+10 

10 

2.973562E+11 

5.983260E-08 

3.342659E+10 

--  Node  »  6  - 

--  Radius  =  5. 

500E-04 - 

1 

1.907419E+17 

2.426457E-08 

0.000000E-f00 

2 

1.152346E+17 

3.471774E-08 

5.760742E+10 

3 

2.970138E+16 

3.980550E-08 

5.024431E+10 

4 

6.723206E+15 

4.386735E-08 

4.559199E+10 

5 

1.396059E+15 

4.730536E-08 

4.227850E+10 

6 

2.723131E+14 

5.031653E-08 

3.974836E+10 

7 

5.061627E+13 

5.301388E-08 

3.772596E+10 

8 

9.051950E+12 

5.546897E-08 

3.605619E+10 

9 

1.577365E+12 

5.773034E-08 

3.464382E+10 

10 

2.973563E+11 

5.983260E-08 

3.342659E+10 

--  Node  -  7  — 

--  Radius  -  6. 

500E-04 - 

1 

1.907419E+17 

2.426457E-08 

0.000000E400 

2 

1.152346E+17 

3.471774E-08 

5.760742E+10 

3 

2.970139E+16 

3.980550E-08 

5.024431E+10 

4 

6.723208E+15 

4.386735E-08 

4.559199E+10 

5 

1.396060E+15 

4.730536E-08 

4.227850E+10 

6 

2.723132E+14 

5.031653E-08 

3.974836E+10 

7 

5.061626E+13 

5.301388E-08 

3.772596E+10 

8 

9.051948E+12 

5.546897E-08 

3.605619E+10 

9 

1.577364E+12 

5.773034E-08 

3.464382E+10 

10 

2.973560E+11 

5.983260E-08 

3.342659E+10 

--  Node  «  8  — 

--  Radius  =  7.000E-04 - 

1 

1.907419E+17 

2.426457E-08 

0.000000E-f«0 

2 

1.152345E+17 

3.471774E-08 

5.760742E+10 

3 

2.970137E+16 

3.980550E-08 

5.024431E+10 

4 

6.723204E+15 

4.386735E-08 

4.559199E+10 

5 

1.396059E+15 

4.730536E-08 

4.227850E+10 

6 

2.723131E+14 

5.031653E-08 

3.974836E+10 

7 

5.061627E+13 

5.301388E-08 

3.772596E+10 

8 

9.051952E+12 

5.546897E-08 

3.605619E+10 

9 

1.577365E+12 

5.773034E-08 

3.464382E+10 

10 

2.973565E+11 

5.983260E-08 

3.342659E+10 

--  Node  -  9  — 

—  Radius  -  7.500E-04 - 

1 

1.907420E+17 

2.426457E-08 

0.000000E-f00 

2 

1.152348E+17 

3.471774E-08 

5.760742E+10 

3 

2.970145E+16 

3.980550E-08 

5.024431E+10 

4 

6.723221E+15 

4.386735E-08 

4.559199E+10 

5 

1.396062E+15 

4.730536E-08 

4.227850E+10 

6 

2.723133E+14 

5.031653E-08 

3.974836E+10 

7 

5.061624E+13 

5.301388E-08 

3.772596E+10 

8 

9.051934E+12 

5.546897E-08 

3.605619E+10 

9 

1.577360E+12 

5.773034E-08 

3.464382E+10 

10 

2.973548E+11 

5.983260E-08 

3.342659E+10 

230 

--  Node  -  10  -  — 

—  Radius  =  8.000E-04 - 

1 

1.907418E+17 

2.426457E-08 

0.000000E-f00 

2 

1.152344E+17 

3.471774E-08 

5.760742E+10 

3 

2.970132E+16 

3.980550E-08 

5.024431E+10 

4 

6.723194E+15 

4.386735E-08 

4.559199E+10 

5 

1.396058E+15 

4.730536E-08 

4.227850E+10 

6 

2.723130E+14 

5.031653E-08 

3.974836E+10 

7 

5.061629E+13 

5.301388E-08 

3.772596E+10 

8 

9.051963E+12 

5.546897E-08 

3.605619E+10 

9 

1.577369E+12 

5.773034E-08 

3.464382E+10 

10 

2.973574E+11 

5.983260E-08 

3.342659E+10 

--  Node  *  11  — 

--  Radius  =  8. 

500E-04 - 

j 

1 

1.907426E+17 

2.426457E-08 

0.000000E+00 

2 

1.152358E+17 

3.471774E-08 

5.760742E+10 

i 

3 

2.970175E+16 

3.980550E-08 

5.024431E+10 

1 

4 

6.723284E+15 

4.386735E-08 

4.559199E+10 

5 

1.396071E+15 

4.730536E-08 

4.227850E+10 

6 

2.723141E+14 

5.031653E-08 

3.974836E+10 

7 

5.061614E+13 

5.301388E-08 

3.772596E+10 

8 

9.051868E+12 

5.546897E-08 

3.605619E+10 

9 

1.577339E+12 

5.773034E-08 

3.464382E+10 

10 

2.973490E+11 

5.983260E-08 

3.342659E+10 

--  Node  -  12  — 

--  Radius  -  8. 

900E-04 . . 

! 

1 

1.907413E+17 

2.426457E-08 

0.0000e0E400 

j 

2 

1.152335E+17 

3.471774E-08 

5.760742E+10 

3 

2.970104E+16 

3.980550E-08 

5.024431E+10 

4 

6.723135E+15 

4.386735E-08 

4.559199E+10 

5 

1.396049E+15 

4.730536E-08 

4.227850E+10 

6 

2.723123E+14 

5.031653E-08 

3.974836E+10 

7 

5.061639E+13 

5.301388E-08 

3.772596E+10 

8 

9.052025E+12 

5.546897E-08 

3.605619E+10 

9 

1.577388E+12 

5.773034E-08 

3.464382E+10 

10 

2.973629E+11 

5.983260E-08 

3.342659E+10 

--  Node  =  13  — 

—  Radius  =  9.300E-04 - 

1 

1.907460E+17 

2.426457E-08 

0.000000E-f00 

2 

1.152416E+17 

3.471774E-08 

5.760742E+10 

3 

2.970360E+16 

3.980550E-08 

5.024431E+10 

4 

6.723668E+15 

4.386735E-08 

4.559199E+10 

5 

1.396128E+15 

4.730536E-08 

4.227850E+10 

6 

2.723188E+14 

5.031653E-08 

3.974836E+10 

7 

5.061556E+13 

5.301388E-08 

3.772596E+10 

8 

9.051475E+12 

5.546897E-08 

3.605619E+10 

9 

1.577216E+12 

5.773034E-08 

3.464382E+10 

10 

2.973139E+11 

5.983260E-08 

3.342659E+10 

—  Node  =  14  — 

--  Radius  =  9. 

500E-04 - 

1 

1.907370E+17 

2.426457E-08 

0.000000E400 

2 

1.152263E+17 

3.471774E-08 

5.760742E+10 

3 

2.969876E+16 

3.980550E-08 

5.024431E+10 

4 

6.722662E+15 

4.386735E-08 

4.559199E+10 

5 

1.395979E+15 

4.730536E-08 

4.227850E+10 

6 

2.723068E+14 

5.031653E-08 

3.974836E+10 

1 
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7  5.061717E+13  5.301388E-08  3.772596E+10 

8  9.052528E+12  5.546897E-08  3.605619E+10 

9  1.577544E+12  5.773034E-08  3.464382E+10 

10  2.974071E+11  5.983260E-08  3.342659E+10 

-  Node  »  15  -  Radius  «  9.700E-04 - 

1  1.907768E+17  2.426457E-08  0.000000E+00 

2  1.152929E+17  3.471774E-08  5.760742E+10 

3  2.971989E+16  3.980550E-08  5.024431E+10 

4  6.727095E+15  4.386735E-08  4.559199E+10 

5  1.396649E+15  4.730536E-08  4.227850E+10 

6  2.723653E+14  5.031653E-08  3.974836E+10 

7  5.061180E+13  5.301388E-08  3.772596E+10 

8  9.048326E+12  5.546897E-08  3.605619E+10 

9  1.576200E+12  5.773034E-08  3.464382E+10 

10  2.970200E+11  5.983260E-08  3.342659E+10 

- Node  «  16 - Radius  =  9.800E-04 - 

1  1.907013E+17  2.426457E-08  0.000000E+00 

2  1.151657E+17  3.471774E-08  5.760742E+10 

3  ;.967965E+16  3.980550E-08  5.024431E+10 

4  6.718704E+15  4.386735E-08  4.559199E+10 

5  1.395400E+15  4.730536E-08  4.227850E+10 

6  2.722627E+14  5.031653E-08  3.974836E+10 

7  5.062442E+13  5.301388E-08  3.772596E+10 

8  9.056896E+12  5.546897E-08  3.605619E+10 

9  1.578883E+12  5.773034E-08  3.464382E+10 

10  2.977862E+11  5.983260E-08  3.342659E+10 

. Node  -  17 . Radius  »  9.900E-04 . 

1  1.910678E+17  2.426457E-08  0.000000E+00 

2  1.157624E+17  3.471774E-08  5.760742E+10 

3  2.987072E+16  3.980550E-08  5.024431E+10 

4  6.759842E+15  4.386735E-08  4.559199E+10 

5  1.401990E+15  4.730536E-08  4.227850E+10 

6  2.729704E+14  5,031653E-08  3.974836E+10 

7  5.062448E+13  5.301388E-08  3.772596E+10 

8  9.030537E+12  5.546897E-08  3.605619E+10 

9  1.569303E+12  5.77303-1E-08  3.464382E+10 

10  2.948855E+11  5.983260E-08  3.342659E+10 

-  Node  -  18  -  Radius  =  9.925E-04 - 

1  1.912319E+17  2.426457E-08  0.000000E+00 

2  1.154131E+17  3.471774E-08  5.760742E+10 

3  2.983335E+16  3.980550E-08  5.024431E+10 

4  6.792335E+15  4.386735E-08  4.559199E+10 

5  1.421386E+15  4.730536E-08  4.227850E+10 

6  2.797620E+14  5.031653E-08  3.974836E+10 

7  5.251273E+13  5.301388E-08  3.772596E+10 

8  9.488197E+12  5.546897E-08  3.605619E+10 

9  1.671186E+12  5.773034E-08  3.464382E+10 

10  3.188721E+11  5.983260E-08  3.342659E+10 

-  Node  -  19  -  Radius  -  9.950E-04 - 


1  2.012295E+17  2.426457E-08  0.000000E400 
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2  1.291832E+17  3.471774E-08  5.760742E+10 

3  3.460673E+16  3.980550E-08  5.024431E+10 

4  8.009758E+15  4.386735E-08  4.559199E+10 

5  1.682402E+15  4.730536E-08  4.227850E+10 

6  3.296526E+14  5.031653E-08  3.974836E+10 

7  6.127661E+13  5.301388E-08  3.772596E+10 

8  1.092718E+13  5.546897E-08  3.605619E+10 

9  1.895029E+12  5.773034E-08  3.464382E+10 

10  3.547876E+11  5.983260E-08  3.342659E+10 

-  Node  =  20  -  Radius  -  9.987E-04 - 


1  1.107092E+17  2.426457E-08  0.000000E+00 

2  4.447289E+16  3.471774E-08  5.760742E+10 

3  7.309087E+15  3.980550E-08  5.024431E+10 

4  1.026271E+15  4.386735E-08  4.559199E+10 

5  1.281519E+14  4.730536E-08  4.2278S0E+10 

6  1.434489E+13  S.031653E-08  3.974836E+10 

7  1.403750E+12  S.301388E-08  3.772596E+10 

8  1.078410E+11  5.546897E-08  3.605619E+10 

9  3.013661E+09  5.773034E-08  3.464382E+10 

10  0.000000E-f00  5.983260E-08  3 . 342659E-t-10 

-  Node  =  21  -  Radius  =  1.000E-03 - 

1  0.000000E4«0  2.42e457E-08  0.000000E-^ 

2  0.000000E+00  3.471774E-08  5.760742E+10 

3  0.000000E^  3.980550E-08  5.024431E-h10 

4  0.000000E+00  4.386735E-08  4.559199E+10 

5  0.000000E+00  4.730536E-08  4 . 227850E-t-10 

6  0.000000E400  5.031653E-08  3.974836E+10 

7  0.000000E+00  5.301388E-08  3.772596E+10 

8  0.000000E-f00  5.546897E-08  3.605619E-fl0 

9  0.000000E+00  5.773034E-08  3.464382E+10 

10  0.000000E+00  5.983260E-08  3.342659E+10 
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